Fast hierarchical tomography methods and apparatus

ABSTRACT

Pixel images ƒ are created from projections (q 1  . . . q p ) by backprojecting ( 100 ) selected projections to produce intermediate images (I 1 , m), and performing digital image coordinate transformations ( 102 ) and/or resampling (FIG.  31, 186, 192, 196 ) on selected intermediate images. The digital image coordinate transformations ( 102 ) are chosen to account for view angles of the constituent projections of the intermediate images and for their Fourier characteristics, so that the intermediate images may be accurately represented by sparse samples. The resulting intermediate images are aggregated into subsets ( 104 ), and this process is repeated in a recursive manner until sufficient projections and intermediate images have been processed and aggregated to form the pixel image ƒ. Digital image coordinate transformation can include rotation (FIG.  18, 102 ), shearing (Fig. IOB,  120, 122 ), stretching, contractions ( 109 ), etc. Resampling can include up-sampling ( 101, 106 ), down-sampling ( 109 ), and the like. Projections (FIG.  32,  pθ 1  . . . pθ 18 ) can be created from a pixel image (f), by performing digital iamge coordinate transformation ( 202 ) and/or resampling ( 204 ) and/or decimation (FIG.  32, 204;  FIG.  33, 212 ) re-projecting the final intermediate image ( 208 ).

This is a continuation-in-part of Provisional Application Ser. No.60/501,350, filed Sep. 9, 2003, incorporated by reference in itsentirety.

This invention relates to tomography, and more particularly, to methodsand apparatus for creating pixel images from projections.

BACKGROUND OF THE INVENTION

Tomographic reconstruction is a well-known technique underlying nearlyall of the diagnostic imaging modalities including x-ray computedtomography (CT), positron emission tomography (PET), singly photonemission tomography (SPECT), and certain acquisition methods formagnetic resonance imaging (MRI). It also finds application inmanufacturing for nondestructive evaluation (NDE), for securityscanning, in synthetic aperture radar (SAR), radio astronomy, geophysicsand other areas.

There are several main formats of tomographic data: (i) parallel-beam,in which the line-integrals are performed along sets of parallel lines;(ii) divergent-beam, in which the line-integrals are performed alongsets of lines that diverge as a fan or a cone; and(iii) curved, in whichthe integrals are performed along sets of curves, such as circles,ellipses, or other closed or open curves. One problem of tomographicreconstruction is to reconstruct a 2D or 3D image from a set of itsline-integral projections. Another problem of tomographic reconstructionis to reconstruct a 3D image from a set of its surface-integralprojections, that is, its integrals on a family of surfaces. Forexample, the 3D Radon transform involves integrals of the image on afamily of 2D planes of various orientations and distances from theorigin. Some of the problems of tomographic reconstruction, and some ofthe reconstruction methods, are described in standard references such asF. Natterer, The Mathematics of Computerized Tomography. Chichester:John Wiley, 1986; F. Natterer and F. Wubbeling, Mathematical Methods inImage Reconstruction. Philadelphia: Society for Industrial and AppliedMathematics, 2001; A. C. Kak and M. Soaney, Principles of ComputerizedTomographic Imaging. New York: IEEE Press, 1988; and S. R. Deans, TheRadon Transform and Some of Its Applications. New York: Wiley, 1983.

The method of choice for tomographic reconstruction is filteredbackprojection (FBP) or convolution backprojection (CBP), which use anunweighted (in the parallel-beam or Radon Transform cases)or a weighted(in most other cases) backprojection step. This step is thecomputational bottleneck in the technique, with computationalrequirements scaling as N³ for an N×N-pixel image in 2D, and at least asN⁴ for an N×N×N-voxel image in 3D. Thus, doubling the image resolutionfrom N to 2N results in roughly an 8-fold (or 16-fold, in 3D) increasein computation. While computers have become much faster, with the adventof new technologies capable of collecting ever larger quantities of datain real time (e.g., cardiac imaging with multi-row detectors,interventional imaging), and the proliferation of 3D acquisitiongeometries, there is a growing need for fast reconstruction techniques.Fast reconstruction can either speed up the image formation process,reduce the cost of a special-purpose image reconstruction computer, orboth.

The dual operation of backprojection is reprojection, which is theprocess of computing the projections of an electronically stored image.This process, too, plays a fundamental role in tomographicreconstruction. A combination of backprojection and reprojection canalso be used to construct fast reconstruction algorithms for the longobject problem in the helical cone-beam geometry, which is key topractical 3D imaging of human subjects. Furthermore, in variousapplications it is advantageous or even necessary to use iterativereconstruction algorithms, in which both backprojection and reprojectionsteps are performed several times for the reconstruction of a singleimage. Speeding up the backprojection and reprojection steps willdetermine the economic feasibility of such iterative methods.

Several methods have been proposed over the years to speed upreconstruction. For example, Brandt et al. U.S. Pat. No. 5,778,038describes a method for 2D parallel-beam tomography using a multileveldecomposition, producing at each stage an image covering the entirefield-of-view, with increasing resolution. Nillson et al. U.S. Pat. No.6,151,377 disclose other hierarchical backprojection methods. Whilethese systems may have merit, there is still a need for methods andapparatus that produce more accurate images, and offer more flexibilitybetween accuracy and speed.

Accordingly, one object of this invention is to provide new and improvedmethods and apparatus for computed tomography (CT) scanning.

Another object is to provide new and improved methods and apparatus forCT scanning that produce more accurate images, and offer moreflexibility between accuracy and speed.

SUMMARY OF THE INVENTION

These objects are achieved or exceeded by the present invention. Pixelimages are created from projections by backprojecting selectedprojections to produce intermediate images, and performing digital imagecoordinate transformations and/or resampling on intermediate images. Thedigital image coordinate transformations are chosen to account for viewangles of the constituent projections of the intermediate images and fortheir Fourier characteristics, so that the intermediate images may beaccurately represented by sparse samples. The resulting intermediateimages are aggregated into subsets, and this process is repeated in arecursive manner until sufficient projections and intermediate imageshave been processed and aggregated to form the pixel image.

Digital image coordinate transformation can include rotation, shearing,stretching, contractions, etc. Reampling can include up-sampling,down-sampling and the like.

Projections can be created from a pixel image by performing digitalimage coordinate transformation and/or resampling and/or decimation andre-projecting the final intermediate image.

BRIEF DESCRIPTION OF THE DRAWINGS

The above mentioned and other features of this invention and the mannerof obtaining them will become more apparent, and the invention itselfwill be best understood by reference to the following description of anembodiment of the invention taken in conjunction with the accompanyingdrawings, in which:

FIG. 1 is a block diagram of apparatus used for the present invention;

FIGS. 2A, 2B and 2C are diagrams of sampling patterns used in someembodiments of the present invention;

FIGS. 3A, 3B and 3C are additional sampling patterns used in someembodiments of the present invention;

FIG. 4 is a diagram illustrating a known method of backprojection;

FIG. 5A is a diagram illustrating an algorithm for one embodiment of thepresent invention;

FIG. 5B is a diagram illustrating the manner in which intermediateimages are generated in the embodiment of FIG. 5A;

FIG. 6 is a diagram illustrating Fourier characteristics used to producethe intermediate images of FIG. 5A;

FIGS. 7A, 7B and 7C are diagrams showing Fourier supports ofintermediate images for the backprojection algorithm illustrated in FIG.5A, when the coordinate transformation is a digital image rotation;

FIG. 8 is a diagram illustrating an algorithm used in another embodimentof the present invention;

FIG. 9 is a diagram showing the evolution of the spectral support in thealgorithm of FIG. 8, the blocks (1 . . . 9) corresponding to thecorresponding points in FIG. 8;

FIG. 10A is a diagram describing an algorithm for embodiment of thepresent invention;

FIG. 10B illustrates the image coordinate transformation used in theembodiment of FIG. 10A;

FIG. 11A is a diagram illustrating shear scale backprojection, and FIG.11B is a diagram illustrating hierarchical shear scale backprojection;

FIGS. 12A and 12B are diagrams showing the effect of image shearing onthe spectral support of intermediate images;

FIG. 13 is a diagram illustrating an algorithm for another embodiment ofthe present invention;

FIG. 14 is an algorithm for finding optimal shear factors;

FIG. 15 illustrates an algorithm for another embodiment of the presentinvention;

FIG. 16 illustrates an algorithm for still another embodiment of thepresent invention;

FIG. 17 is a diagram which illustrates common fan beam geometry with acircular scanning trajectory;

FIG. 18 illustrates an algorithm for another embodiment of the presentinvention;

FIG. 19 illustrates weighting of intermediate images in the algorithmdescribed in FIG. 18;

FIGS. 20A, 20B and 20C illustrate sampling points for the secondhierarchical level of FIG. 18;

FIGS. 21A and 21B are diagrams of sampling patterns used in thealgorithm of FIG. 18;

FIG. 22 is a diagram showing original intersection points obtained usingthe method illustrated in FIGS. 20A-20C;

FIGS. 23A, 23B, 23C and 23D illustrate sampling points for rotation forup-sampling used in FIG. 18;

FIG. 24 illustrates local spectral support at a point of an intermediateimage generated by the algorithm of FIG. 18;

FIG. 25 is a diagram of nonuniform sampling patterns used in thealgorithm of FIG. 18;

FIGS. 26A and 26B illustrate sampling patterns for resampling andcoordinate transformation;

FIGS. 27A and 27B illustrate alternative sampling schemes which can beused in the present invention;

FIG. 28 includes two diagrams of divergent beams used in the presentinvention;

FIG. 29A is a diagram showing a conebeam, and FIG. 29B is a diagramillustrating resampling;

FIG. 30 is a diagram of an algorithm used for resampling projections;

FIG. 31 is another algorithm used for resampling in the embodiment ofthe present invention;

FIG. 32 is a diagram of diagram of an algorithm used for fasthierarchical reprojection;

FIG. 33 is a diagram of another algorithm for fast hierarchicalreprojection;

FIG. 34 is a graph showing the results of experiments using the presentinvention;

FIG. 35 includes sample images generated with the present invention;

FIG. 36A is a display of a reconstructed image using the conventionalalgorithm, and FIG. 36B shows a result obtained with the fast algorithmsof the present invention; and

FIGS. 37A and 37B are diagrams of point spread functions of algorithmscomparable to conventional algorithms, and FIG. 37C displays a pointspread function of a fast algorithm of the present invention.

DETAILED DESCRIPTION Symbols and Fonts

The following system of mathematical symbols and fonts will be used toimprove clarity.

Functions in the space domain are denoted by small letters (e.g.n ƒ(x)), while their Fourier transforms are denoted by capital letters(F(ω)).

The indices of two-variable functions are denoted variously, dependingon convenience. The following three notations of function ƒ areequivalent: ƒ( x₁, x₂), ƒ( {right arrow over (x)}), and${f\left( \begin{bmatrix}{x\quad 1} \\{x\quad 2}\end{bmatrix} \right)}.$

Continuous-domain and discrete-domain functions respectively aredistinguished by the style of parentheses used with their indices: ƒ(x₁, x₂) is a function of two continuous variables (i.e., {right arrowover (x)} ε R²), and f[m₁, m₂] is the sampled version of ƒ( {right arrowover (x)}) and is therefore a 2-D array.

A linear operator and its corresponding matrix are distinguished by thefont style. Suppose (A ε R²x²) is a matrix, then A is its associatedlinear operator. For example, if A is a coordinate transformation,g({right arrow over (x)})=(Af)({right arrow over (x)})=ƒ( A{right arrowover (x)}).

The same operator is sometimes denoted differently inside and outsideblock diagrams. While outside it may be denoted as A(α), within theblock diagram it is denoted as A_(α).

Overview of Hardware

The present invention has application in a variety of imaging apparatus,including CT scanners. Typical imaging apparatus 10 (FIG. 11) includes ascanner 12 which acquires data from an object such as a head, and sendsraw data corresponding to line-integral projections, e.g., with adivergent beam geometry, 14 to a projection pre-processor 16. Theprojection pre-processor 16 applies various conversions, normalizations,and corrections to the data, as well as weighting and filtering, whichmay be shift varying. The output of the projection pre-processor 16 is acollection of pre-processed projections, hereinafter simply referred toas projections, also called sinogram, 18, which is fed to a sinogramupdate processor 20. The sinogram update processor 20 possibly modifiesthe input sinogram₁ 18, using information from sinogram₂ 34, for examplecorrecting for various artifacts including beam-hardening, or as part ofa multi-step or iterative reconstruction procedure.

The output of the sinogram update processor 20 is a sinogram₃ 22, whichis input to a fast backprojection processor 24. The fast backprojectionprocessor 24 is generally a computer or special purpose hardware, or anycombination thereof, of any suitable type programmed and/or wired toperform the algorithms described herein.

The output of the fast backprojection processor 24 is an electronicimage, 26, which is input to an image conditioning processor 28. Theimage conditioning processor 28 performs necessary postprocessing of theelectronic image, possibly including identification and extraction ofartifact images, or images for further processing in a multi-step oriterative reconstruction process.

If desired, the image conditioning processor 28 can produce anelectronic image₂ 30 that is fed to a fast reprojection processor 32.The fast reprojection processor 32 is generally a computer or specialpurpose hardware, or any combination thereof, of any suitable typeprogrammed and/or wired to perform the algorithms described herein. Ifdesired, this processor can share the use of the same computer andhardware employed by the backprojection processor 24.

The output of the fast reprojection processor 32 is a sinogram₂ 34,which is fed back into the sinogram update processor 20. Thebackprojection/reprojection process can continue until suitable resultsare obtained. While reprojection is not always needed, it is helpful inmany situations.

When the electronic image, 26 is suitable, the image conditioningprocessor 28 produces an electronic image₃ 36, which is fed to astorage/analysis/display device 38. It is contemplated that theelectronic image₃ 36 can be stored in a computer memory, and/or analyzedelectronically for anomalies or dangerous materials, for example, and/ordisplayed, and/or printed in some viewable form.

Overview of Backprojection and Reprojection Methods of the PresentInvention

The backprojection methods of the present invention use varioustechniques to create an image made of pixels (picture elements) and/orvoxels (3D picture elements), hereinafter referred to collectively aspixels, which will now be introduced in a general way.

This explanation uses terminology and processes commonly used inmulti-dimensional signal processing, for example as described in D.Dudgeon and R. Mersereau, Multidimensional Digital Signal Processing.Englewood Cliffs: Prentice-Hall, 1983. Some terms in this description ofthe present invention are used in the following contexts. The termSampling pattern refers to a set of points in space with positionsdefined relative to a system of coordinates. Examples of samplingpatterns are seen in FIGS. 2A-2C and 3A-3C. A Cartesian sampling patternrefers to a set of points formed by the intersection of two mutuallyperpendicular sets of parallel lines. The term continuous image refersto a function defined on a coordinate system, for example, ƒ( x, y), andƒ( x, y, z) are respectively 2D and 3D functions. A digital image is anarray of values of a continuous image on a sampling pattern. Morebroadly, a continuous image can be represented by an array of numbersthat serve as the coefficients in a series expansion with respect tosome basis set, such as splines, of which the Cartesian product ofzero-th order splines yields the familiar square pixel form fordisplaying digital images as continuous images. Hereinafter, this arrayof numbers will be also referred to as a digital image. All imagesstored in a digital computing device must be digital. For brevity, bothdigital and continuous images will be often referred to hereinaftersimply as images, with the meaning inferred from the context. With thisterminology, a pixel image is a digital image corresponding to asampling pattern that is a lattice, i.e., a uniformly spaced, periodicpattern, usually but not necessarily Cartesian.

One sampling pattern will be said to be sparser than another, if ityields a smaller total number of samples of a continuous image.Typically a sparser sampling pattern will have a lower sampling density.Oversampling refers to using more samples than necessary to represent acontinuous image at a desired accuracy. The corresponding digital imagewill be said to be oversampled. Given a digital image corresponding to acontinuous image for one sampling pattern, the process of producing anew digital image corresponding, with a desired accuracy, to the samecontinuous image on a different sampling pattern, will be called digitalimage resampling. Up-sampling and down-sampling are special cases ofresampling on a denser or sparser sampling pattern, respectively.Further, up-sampling or down-sampling by a factor of 1 involves nochange in the digital image, is considered a form of up-sampling ordown-sampling. Digital image addition refers to to a point-by-pointaddition of digital images defined with respect to the same samplingpattern, or the same basis, in the case of an expansion with respect toa basis. Lower dimensional digital filtering refers to digital filteringof a multidimensional array along only a subset of its dimensions, forexample, separate ID filtering of each column of a 2D rectangulardigital image.

Coordinate transformation of a continuous image refers to operationssuch as rotation, shearing, stretching, or contraction. To define adigital coordinate transformation, consider two continuous imagesrelated by a coordinate transformation, and the corresponding digitalimages representing the continuous images with respect to a commonsampling pattern. The process of producing one digital image from theother is called digital image coordinate transformation. This can beaccomplished by digital filtering, i.e., by discrete index operations onthe digital image. Specific examples include (but are not limited to)digital image rotation, digital image shearing, digital imagestretching, digital image contraction, and the combinations of suchoperations. Methods for performing digital image coordinatetransformation are known, for example, as described in M. Unser, P.Thevenaz, and L. Yaroslavsky, Convolution-based interpolation for fast,high-quality rotation of images. IEEE Transactions Image Processing,Vol. 4, pp. 1371-1381, 1995.

Some digital image coordinate transformations are illustrated in FIGS.2A-2C and 3A-3C. FIG. 2A shows the outline of a continuous image (arectangle) and the sampling pattern for a digital image representing it.Values of the continuous image on the heavy dots are included in thedigital image. FIGS. 2-B and 2-C show the rotated and the shearedcontinuous image, respectively, on the same sampling pattern, with theheavy dots showing the values included in the digitally rotated/shearedversion of the digital image in FIG. 2-A.

FIG. 3A also shows a continuous image and a sampling pattern defining adigital image. FIG. 3B shows the stretched continuous image by someconstant factors in the x and y dimensions. The digitally stretchedimage is defined by values of the stretched continuous image on theheavy dots. FIG. 3C shows the same continuous image as in FIG. 3A, butwith a sampling pattern denser by certain stretch factors. Thecorresponding digital image in FIG. 3C will be the same as in FIG. 3B.More generally, digital image stretching or contraction can beequivalent, for sampling patterns with some regularity, to digital imageup-sampling or down-sampling.

Note that the application of certain coordinate transformations, such asrotation by 0 degrees, shearing by a shear parameter of zero, orstretching or contraction by a factor of 1, leave the digital imageunchanged, and therefore may be either included or omitted in thedescription of a process, without changing the result.

The Fourier transform of a continuous image allows one to determine, viasampling theory, the properties of sampling patterns such that thecorresponding digital image represents the continuous image to desiredaccuracy, as explained, for example, in D. Dudgeon and R. Mersereau,Multidimensional Digital Signal Processing. Englewood Cliffs:Prentice-Hall, 1983. Likewise, the discrete-time Fourier transform(DTFT) of a digital image allows one to determine what digital imagecoordinate transformations will produce a digital image that representsthe transformed continuous image to a desired level of accuracy. Therelevant properties of the Fourier transform of continuous images, andthe DTFT of digital images, will be collectively referred to as Fouriercharacteristics.

Weighted backprojection operations require weighting of each projectionby a weight that depends on the position of the pixel generated by, thebackprojection. Different weights can be used for different projections,depending, for example, on the position of the source at which theprojection was acquired, as explained in A. C. Kak and M. Slaney,Principles of Computerized Tomographic Imaging. New York: IEEE Press,1988. As a special case, the weighting factor can be equal to 1, whichcorresponds to no weighting, but is a weighting factor. Unlessspecifically indicated, weighted and un-weighted backprojection will becollectively referred to as backprojection.

With this background information, several embodiments of the inventionwill be described.

The backprojection processor 28 (FIG. 1) is programmed to perform thealgorithms used to practice the present invention. The algorithms willbe discussed in detail, but will first be described in more generalterms. Steps in the backprojection process are indicated in blockdiagrams, with reference numerals.

FIG. 4 illustrates rotation-based backprojection, as the sum of rotatedintermediate images formed by backprojecting individual projections q₁ .. . q_(p) at zero angle in step 50. The backprojections produces imageswhich are subjected to coordinate transformation at 52, and aggregatedat 54 to produce an image {circumflex over (f)}. By itself, thisstructure is equivalent to conventional backprojection, and offers noreduction in operation count. However, it serves as a stepping stone tothe introduction of some of the fast hierarchical backprojection methodsof the present invention.

FIG. 5A illustrates a hierarchical backprojection method for creating apixel image {circumflex over (f)} from a plurality of projections q₁ . .. q_(p). Each projection q_(m) is backprojected at 100 to produce aplurality of intermediate images I_(1,1) . . . I_(1,p). This is thezeroth or preparatory level of a hierarchical backprojection. Digitalimage coordinate transformations are performed on selected intermediateimages at 102. Subsets of the transformed intermediate images (pairs,here), are aggregated at 104 to produce aggregate intermediate imagesI_(2,1), . . . I_(2,P/2). This is the first level of a hierarchicalbackprojection. The aggregate intermediate images of the first levelserve as new intermediate images input to the next level of hierarchicalbackprojection. The process of applying digital image coordinatetransformations to selected intermediate images at 106, and aggregatingselected intermediate images at 108 to produce new intermediate imagescontinues until all intermediate images have been processed andaggregated into the final image {circumflex over (f)} at 116.

If desired, the operations within and across some of the levels can becombined. For example, the zeroth and first level can be combined forsome of the projections, and the corresponding intermediate images fromthe set I_(2,1), . . . I_(2,P/2) produced instead by a backprojection112 at selected view angles of two or more (exactly two, for theembodiment in FIG. 5A) selected projections q_(p), as shown in FIG. 5B.Alternatively, some of the initial intermediate images can be producedby an equivalent process involving no explicit backprojection, such as adirect Fourier reconstruction method, as described in F. Natterer and F.Wubbeling, Mathematical Methods and Image Reconstruction. Philadelphia:Society for Industrial and Applied Mathematics, 2001.

The parameters of the various digital image coordinate transformationsare chosen to account for the view angles of the constituent projectionsof the intermediate images, and for the Fourier characteristics of theintermediate images, so that the aggregates of the intermediate imagesmay be accurately represented by sparse samples, as will be explained.These Fourier characteristics focus on the essential spectral support ofthe intermediate images, i.e. the region in the Fourier (frequency)domain, where the Fourier transform of the intermediate image issignificantly different from zero. Sampling theory teaches that thespectral support of a continuous image determines the nature of thesampling patterns that produce digital images that represent theunderlying continuous image, and from which the continuous image can bereliably reconstructed. In particular, FIG. 6 shows the typicalwedge-shaped spectral support of intermediate image Ĩ_(l,m) in thehierarchical algorithm.

FIGS. 7(A) and 7(B) show Fourier supports of intermediate images in thebinary hierarchical backprojection algorithm illustrated in FIG. 5(A),when the coordinate transformation is chosen to be a digital imagerotation. FIG. 7(A) shows the Fourier-domain support of the virtualintermediate images Ĩ. A virtual image Ĩ_(i,m) is composed of thebackprojection of the projections that are included in the correspondingimage I_(l,m). FIG. 7(B) shows that by choosing the parameters of thecoordinate transformation to account for the view angles of theconstituent projections of the intermediate images, and for the Fouriercharacteristics of the intermediate images, the vertical bandwidth (theheight of the broken-line rectangle) of the intermediate images can beminimized, allowing sparse sampling and reducing computationalrequirements. FIG. 7(C) shows the space-domain sampling scheme ofI_(l,m),I_(l−1,2m) and I_(l−1,2m−1). The sampling points are at theintersections of the dotted lines in the space domain.

FIG. 8 describes another embodiment of the present invention, whichperforms a ternary aggregation of intermediate images, with P=2×3^(L)projections. For concreteness of illustration, the case of P=18projections is shown. As in FIG. 5(A), projections q_(θ) ₁ , q_(θ) ₂ , .. . q_(θ) ₁₈ are backprojected at 100 to produce a plurality ofintermediate images I_(1,1) ^(d) . . . I_(1,18) ^(d). This is the zerothlevel of a hierarchical backprojection.

The projections are grouped in 3's, and selected projections (two of thethree in FIG. 8) in each group are subjected to digital image coordinatetransformation at 102. Subsets (triplets, in FIG. 8) of the transformedintermediate images I_(l) ^(d), are aggregated at 104 to produceaggregate intermediate images I_(2,1) ^(d), . . . , I_(2,6) ^(d). Thisis the first level of hierarchical backprojection (l=1).

In the second level of hierarchical backprojection (l=2) selectedaggregate intermediate images I_(2,m) ^(d) undergo a coordinatetransformation composed of stretching or upsampling along the ycoordinate at 106, and another coordinate transformation K_(2,m) onselected aggregate intermediate images at 108. Here, also, thetransformed intermediate images are aggregated in groups of three at 110to produce intermediate images T_(3,1) ^(d), I_(3,2) ^(d). Selectedintermediate images from level two are subjected to digital imagecoordinate transformations at 112 and aggregated at 114 to produce thenext level intermediate images I_(4,1) ^(d). This process is repeated tothe extent necessary to produce the image {circumflex over (f)} at 116,depending on the number of projections. The digital image coordinatetransformation denoted by K at the last level 112 is a digital imagerotation, and those at the preceding levels 102 and 108 can also bechosen to be digital image rotations. Here too, the parameters of thevarious digital image coordinate transformation are chosen to accountfor the view angles of the constituent projections of the intermediateimages, and for the Fourier characteristics of the intermediate images,so that the aggregates of the intermediate images may be accuratelyrepresented by sparse samples, as will be explained.

FIG. 9 shows the evolution of the spectral support in the ternarybackprojection algorithm of FIG. 8. The numbers (1), . . . , (9)correspond to the corresponding points in the block diagram of thealgorithm in FIG. 8. The spectra shown are the discrete-time Fouriertransforms (DTFTs) of the digital images I_(l,m).

FIG. 10(A) describes an algorithm for another embodiment of the presentinvention, using ternary two-shear-based hierarchical backprojection.The embodiment of FIG. 10(A) is similar to the embodiment of FIG. 8, andthe reference numerals from FIG. 8 are used where appropriate. As inFIG. 8, P=2×3^(L) projections, and the case of L=2 is shown. Theembodiment of FIG. 10(a) differs from that described in FIG. 8 in tworespects. First, an additional upsampling step along the x coordinate isincluded at 101 in the first level coordinate transformations. Second,at least some of the digital image coordinate transformations, denotedby K at 102 and 108, are composed of a sequence of two digital imageshear operations, the first (120) along the y coordinate (122), thesecond along the x coordinate, as shown in FIG. 10(b).

FIG. 11A illustrates shear-scale backprojection, and FIG. 11Billustrates an equivalent hierarchical shear-scale backprojection. InFIG. 11A, the plurality of projections q₁, . . . , q₄, are backprojectedat 140 and subjected to a shear-scale coordinate transformation at 142.The resulting intermediate images are aggregated at 144.

In FIG. 11B, the projections q₁, . . . , q₄ are backprojected at 146,subjected to a shear-scale coordinate transformation at 148, andaggregated in subsets at 150. The intermediate images produced by theaggregation are subjected to a shear-scale transformation again at 152,and the resulting intermediate images are aggregated at 154. Thisprocess continues in a recursive manner until all of the projections andintermediate images have been processed and aggregated to form an image{circumflex over (f)}. Here too, the parameters of the digital sheartransformations and sampling of intermediate images are selected toaccount for the view angles and the Fourier characteristics, so as toreduce the sampling requirements and required computation.

FIG. 12 shows the effect of image shearing on the spectral support ofintermediate images. FIG. 12A shows the spectral support of a certainsubset of the projections as they should appear in the final image. FIG.12(B) shows the spectral support of the intermediate image composed ofthe same projections, with coordinate transformation parameters chosento minimize the highest radian frequency in the vertical direction.

FIG. 13 illustrates an algorithm for another embodiment of the presentinvention, a ternary hierarchical shear-scale backprojection.Projections q₇₄ ₁, . . . , q_(θ) ₁ , are backprojected at 160 andsubjected to a shear-scale digital image coordinate transformation at162. Subsets of the resulting images are aggregated at 164, and theresulting intermediate images are subjected to up-sampling at 166 and168 digital shear coordinate transformations at 168. Subsets of thoseimages are aggregated at 170, and selected resulting images aresubjected to additional coordinate transformation at 172. This processcontinues until all of the projections and intermediate images have beenprocessed and aggregated at 174, to produce the image {circumflex over(f)}. The final coordinate transformation K_(−π/2) shown here at 172,only involves rearrangement of pixels, when the sampling pattern for thepixel image {circumflex over (f)} is Cartesian.

FIG. 14 is an algorithm for finding parameters of coordinatetransformations for previously described embodiments, using theproperties of the Fourier properties of intermediate images shown inFIG. 12.

FIG. 15 illustrates an algorithm for another embodiment of the presentinvention, over-sampled ternary two-shear hierarchical backprojection.The embodiment of FIG. 15 is similar to the embodiment of FIG. 10A, andthe reference numerals from FIG. 10A are used where appropriate. Inaddition to the steps of FIG. 10A, however, the embodiment of FIG. 15includes a downsampling step 109 in the one before the last, which isthe second level in FIG. 15. Here too, the parameters of the variousdigital image coordinate transformation are chosen to account for theview angles of the constituent projections of the intermediate images,and for the Fourier characteristics of the intermediate images, so thatthe aggregates of the intermediate images may be accurately representedby sparse samples. However, certain degrees of oversampling are used toimprove the accuracy of subsequent processing, as will be explained. Ifdesired, for improved computational efficiency the downsampling step 109can be combined with the second, x-coordinate shear comprising thetwo-shear digital image coordinate transformation 108 (shown in FIG.10B) producing a shear-scale transformation, so that the processes 108and 109 are together replaced by the process shown in FIG. 15(B).

FIG. 16 illustrates an algorithm for another embodiment of the presentinvention, over-sampled ternary hierarchical shear-scale backprojection.The embodiment of FIG. 16 is similar to the embodiment of FIG. 13, andreference numerals from FIG. 13 are used where appropriate. However,FIG. 16 includes additional steps of upsampling 161 in the intermediatelevels, and downsampling 169 in the level before the last level, inwhich the intermediate images are upsampled and downsampled,respectively. Similarly to the embodiment of FIG. 15, in the embodimentof FIG. 16 when downsampling step 169 along the x coordinate ↓_(x)U_(l,m) follows an x-shear step 168 S_(x) ^(α) ^(l,m) , the two can becombined for computational efficiency into a single digital imageshear-scale.

Non-Cartesian Sampling Schemes for Fast Hierarchical BackprojectionAlgorithms

The intermediate images in the different embodiments of hierarchicalbackprojection illustrated in FIGS. 8,10, and 13 have a peculiarspectral support amenable to efficient non-cartesian sampling. Inparticular, the underlying continuous domain image I_(l,m) in the l^(th)level occupies a wedge in Fourier space, as seen for example in FIGS. 6,7, 9, and 12. Multi-dimensional sampling theory tells us that for imageswith a spectral support such as this, sampling on a cartesian grid isless efficient than sampling on an appropriate non-cartesian grid.Non-cartesian sampling can improve sampling efficiency by packing thecopies of the baseband spectrum more tightly in the Fourier plane. Foran explanation of 2D sampling see [?]. In particular, a quincunxsampling scheme reduces the sizes of intermediate images, and thereforethe computational costs of the algorithm, by a factor of almost 2.Digital image coordinate transformation on periodic non-Cartesiansampling patterns can be executed efficiently using one dimensionalshift-invariant filters. Likewise, all the methods previously describedfor selection of the parameters of digital image transformations applyas well in the case of such sampling patterns. Therefore, all theembodiments previously describe extend to embodiments that use periodicnon-Cartesian sampling patterns.

Variants of the embodiments of the present invention already describedare applicable to 3D backprojection of a variety of forms of 3Dprojections, including the X-ray transform that arises for example inPositron Emission Tomography, and the 3D Radon transform, which arisesin magnetic resonance imaging.

The three-dimensional (3D) X-ray transform is a collection of integralsalong sets of parallel lines, at various orientations, in 3D. Each 3DX-ray projection is a two-dimensional function that can be characterizedby the 3D angle at which the lines are oriented. The block-diagram forhierarchical backprojection for 3D X-ray data is similar to the onespreviously described, such as FIGS. 8, 10, 13, 15 or 16. Theintermediate images in this case are three-dimensional, nottwo-dimensional as in the previously described examples. Eachintermediate image is sampled on a 3D sampling pattern that is sparsestin a direction that is an average of the 3D angles of the constituentprojections. As the algorithm progresses the density of samples alongthis sparse-sampling (slow) direction increases to accommodate theincreasing bandwidth in that direction as explained by the Fourieranalysis of 3D X-ray projections. Consequently at every stage in thealgorithm, before the images are aggregated, each has to be upsampledalong this slow direction. The extra dimension available in the 3Dembodiment also provides more coordinate transformations available foruse in the various steps in the algorithm, such as rotations in 3D. Asin the 2D case, the parameters of these digital image coordinatetransformations can be chosen to account for the constituent view anglesand for the Fourier characteristics of the intermediate images. Thesedigital image coordinate transformations can be decomposed into asequence of one-dimensional operations, such as shears and shear-scales,as previously described. As in the 2D case, oversampling in any subsetof the 3 dimensions may be enforced to improve image quality.

A 3D radon transform projection is a one-dimensional function—acollection of integrals along sets of parallel 2D planes, parameterizedby the displacement of the plane from the origin. The view-angle of theprojection is that of the 3D orientation angle of a vector perpendicularto the set of planes. The block-diagram of the hierarchicalbackprojection of 3D radon projections is as in FIGS. 8, 10, 13, 15 or16. In the first level the projections are Radon backprojected onto the3D image domain. These images are constant along two dimensions, andtherefore need be sampled only on the direction perpendicular to the twoconstant directions. When groups of 3D radon projections are combined,the bandwidth of the aggregate image may increase in one or twodimensions, depending on the view-angles of these constituentprojections. It is therefore necessary to upsample the intermediateimage, possibly in two dimensions, before coordinate transforming it andadding to other intermediate images. As in the 2D case, the coordinatetransformations may be performed separably, may be combined with theupsampling operation, and oversampling may be enforced on theintermediate images.

In the various embodiments of the present invention, digital imagecoordinate transformations and downsampling or upsampling operations maybe performed by a sequence of lower (one) dimensional linear digitalfiltering operations. Furthermore, when the sampling patterns used havesome periodicity, these digital filters can be shift invariant, as willbe described in more detail. For computation efficiency, the digitalfilters can be implemented using recursive filter structures, or usingan FFT, as is known. One way to determine preferred values for thedigital filters is using the theory of spline interpolation, asexplained in M. Unser, A. Aldroubi, and M. Eden, Fast B-splinetransforms for continuous image representation and interpolation, IEEETransactions on Pattern Analysis and Machine Intelligence, Vol. 13, pp.277-285, 1991; M. Unser, A. Aldroubi, and M. Eden, B-spline signalprocessing: Part I—theory, IEEE Transactions Signal Processing, Vol. 41,pp. 821-832, 1993; and M. Unser, A. Aldroubi, and M. Eden, B-splinesignal processing: Part II—efficient design and applications, IEEETransactions Signal Processing, Vol. 41, pp. 834-848, 1993.

The hierarchical backprojection methods of the present invention areapplicable to a wide range of tomographic imaging geometries, such asdivergent beam geometries, including fan-beam and cone-beam, witharbitrary source trajectories. The common fan-beam geometry with acircular scanning trajectory is depicted in FIG. 17. The ray sourcemoves on a circular trajectory of radius D around an image of radius R.A fan-beam projection at source angle β corresponds to line integralmeasurements along a set of rays parametrized by fan angle γ. T_(ST) isthe distance along the ray of the source to the closest edge of theimage-disc and T_(END) is the distance of the source to the farthestedge of the disc.

FIG. 18 illustrates an algorithm for an embodiment of the presentinvention, hierarchical ternary rotation-based backprojection,applicable to fanbeam weighted backprojection. The embodiment of FIG. 18is similar to the embodiment of FIG. 8, and the reference numerals fromFIG. 8 are used where appropriate. The embodiment of FIG. 18 differsfrom that described in FIG. 8 in several respects. First, P=4×3^(L)projections, and the case of L=2 is shown. Second, at the zeroth level,the initial intermediate images I_(l,m) ^(d) are produced from thefanbeam projections by weighted backprojection at 99, denoted here by W.Third, at the last stage, four rather than two intermediate images areaggregated. Fourth, the sampling patterns used in most of the earlylevels of the hierarchy are preferably chosen to be non-Cartesian, aswill be explained. Here too, the last digital image coordinatetransformations only require reordering of image pixels, if a Cartesiansampling pattern is used for the intermediate images I_(3,m) ^(d) in theone-before-last level. Also, as in FIG. 5(B), the operations in thezeroth and first level of the hierarchy can be combined, so that theintermediate images I_(2,m) ^(d) are produced by a direct weightedfanbeam backprojection of their three constituent projections, or byother means.

In a selected number of levels, it is beneficial to modify theembodiment of FIG. 18, by including additional weighting steps beforeand after the cascade of upsampling step and digital image rotation(steps 106 and 108 in FIG. 18), as shown in FIG. 19. The intermediateimage I_(2,m) is weighted by spatially varying weights at 180 and 182,respectively before and after steps 106 and 108. As will be explained,this weighting can be used to reduce the sampling requirements of theintermediate images, and thus reduce the computation.

The sampling scheme of the intermediate images afects the performance ofthe algorithm. The desired sampling scheme is one that uses the fewestsamples without losing image information. Here too, the parameters ofthe various digital image coordinate transformations and resamplingoperations can be chosen to account for the view angles of theconstituent projections of the intermediate images, and for the Fouriercharacteristics of the intermediate images, so that the aggregates ofthe intermediate images may be accurately represented by sparse samples,as will be explained. An alternative method for choosing theseparameters is based on intersection of particular sets of rays orcurves, as will be described.

FIGS. 20A-20C illustrate the progression of intermediate images throughthe levels of the ternary hierarchical algorithm described in FIG. 18,for the case of source angles uniformly spaced at intervals Δ_(β). Thefans of the projections that make up intermediate images in thealgorithms are shown. At the zeroth level, each intermediate imageT_(1,m) ^(d) is made up of a single projection with fan oriented at β=0,shown in FIG. 20(A). After the first level of the recursion, eachintermediate image is made up of three projections, with fans as shownin FIG. 20(B) After the second level, each image is made up of ninebackprojected fans, as shown in FIG. 20(C).

The intersection-based method for choosing coordinate transformationsand sampling/resampling patterns in the algorithm is illustrated inFIGS. 21A amd 21B, for intermediate image I_(3,m), with constituent fansas shown in FIG. 20(C). The sampling pattern for I_(3,m) is made up ofpoints that lie on the central constituent fan, which coincides with theone shown in FIG. 20(A). As shown in FIGS. 21(A) and 21(B), the samplingpoints for I_(3,m) are determined by the intersection of half of thecentral fan with the extremal constituent fans on the respective side ofthe central fan.

It is advantageous to modify the resulting sampling pattern by applyingtwo constraints, to improve the accuracy of the backprojection, andreduce the computation requirements. First, the density of samples alongthe rays is limited not to exceed the desired sampling density of thefinal image. Second, the sampling pattern is forced to contain on eachray at least one sample on the outside of the image disk on each side.The plots in FIG. 22 displays the position of sampling points along aparticular ray of the central fan at fan angle γ_(r). The samplingpoints lie at constant intervals in γ′, where γ′ is the fan angle of oneof the extramal fans shown in FIGS. 21(a)(b) and (c). The points thatfall on the continuous curve in FIG. 22 are original intersection pointsobtained using the method illustrated in FIGS. 21A and 21B. The pointson the broken curve are those modified using the above two constraints.?FIG. 23(D)? shows an example of a sampling pattern obtained using thismodified intersection method. The fan shown is the central fan of theintermediate image for which this sampling pattern has been produced.

As in the case of the previously described embodiment of the presentinvention, it is advantageous to decompose the resampling and digitalimage rotations operations into a sequence of one dimensionaloperations. The blocks marked ↑ U in FIG. 18 represent the upsampling ofthe intermediate image sampled on a central fan onto a finer set ofsampling points on the same fan. This can be achieved by separateupsampling along each ray of the fan. For the cascade of upsampling androtation marked by ↑ U and K in FIG. 18, it is conveniet to do thedecomposition into 1D operation jointly, as illustrated in FIGS.23(A)-23(D). In each of the four panels, the two dashed circlesrepresent the boundary of the image of radius D and the sourcetrajectory on the (larger)circle of radius R. The sampling points arerepresented by small circles. The digital intermediate image withcentral fan (Fan_(a)) and the sampling pattern shown in FIG. 23(A) is tobe resampled and rotated to the angle of the central fan (Fan_(d)) shownin FIG. 23(D), with the final sampling pattern shown in FIG. 23(d). Thisis accomplished in the following two steps:

-   -   (i) Upsample the image, separately along each ray of Fan_(a), to        the sampling pattern shown in 23(b), which is defined by the        intersection of Fan_(a) with Fan_(d). These points will        therefore lie on Fan_(d), as shown in FIG. 23(c);    -   (ii) Upsample the image separately along each ray of Fan_(d), to        the sampling pattern shown in FIG. 23(d).

These steps accomplish the combined operations of resampling androtation, using 1D upsampling operations.

The Fourier-analysis techniques described for the embodiments of thepresent invention in the case of intermediate images sampled on periodicsampling patterns are extended to devise spatially varying samplingschemes in the divergent beam case. These techniques are general enoughto be applied to projections and back-projections on arbitrarytrajectories, in both two and three dimensions. For the case ofbackprojection on non-periodic systems of lines or curves, the conceptof local spectral support replaces that of spectral support. This isillustrated in FIG. 24 for an intermediate image produced by thebackprojection of a single fan shown on the left side of FIG. 24. Aswill be explained, the local spectral of this continuous intermediateimage at the indicated point parametrized by r and θ is the line segmentin the Fourier domain shown on the right in FIG. 24(b). The localspectral support at a point of an intermediate image composed of thebackprojection of multiple fanbeam projections is shown in FIG. 25. Onthe left, the position of the point is indicated, as well as the rangeof view angles of the constituent projections for the intermediateimage. The bow-tie shaped region on the right is the corresponding localspectral support.

This analysis of the local spectral support is used to determine localsampling requirements for the intermediate image. The resultingspatially nonuniform sampling patterns are indicated by the dotted arcsin FIGS. 26A and 26B for two typical intermediate fan-backprojectedimages. The image boundary is indicated by the circle in broken line,and only few points need be taken outside this boundary. This localFourier sampling method may be applied directly to find sampling schemesfor arbitrary projection geometries over lines, curves or planes overarbitrary dimensions.

The Fourier-based method for determining the sampling patterns forresampling and coordinate transformation for hierarchical fan-beambackprojection can be further extended, as will be explained, tosampling patterns lying on lines or curves other than central fans ofthe intermediate images. Examples of such beneficial sets areillustrated in FIGS. 27(A) and 27(B).

General Divergent-Beam Algorithms

The methods described for fanbeam backprojection, extend directly toother divergent-beam geometries, including 3D cone-beam, which is one ofthe most important in modern diagnostic radiology as will be described.Similarly to the fan-beam geometry, in which a source of divergent raysmoves on a trajectory (i.e., a circle) around the imaged object, in thegeneral divergent-beam geometry a source of divergent rays moves on a(not necessarily circular) trajectory in 3D space around the imagedobject. The detector surface measures line integrals through the imagedobject.

One embodiment of a hierarchical backprojection algorithm for a generaldivergent-beam geometry can again be described by a block diagramsimilar to FIG. 18, but modified in the following ways. First, at thezeroth level, the divergent-beam projections are zero-backprojected at99 with the appropriate divergent-beam single-view backprojection W_(l)corresponding to a nominal “zero” source position, producing initialintermediate images. Second, because the trajectory of the source is notnecessarily circular, the constituent divergent-beams of theintermediate images may not simply rotated as in the fan-beam geometry,but also translated, with respect to each other. The coordinatetransformations K in 18 are selected accordingly. Third, depending onthe presence of symmetries in the source trajectory and positions, theremay or may not be “free” coordinate transformations such as the pixelre-arrangement which replaces a π/2 rotation in the fan-beam algorithm.

As in the fan-beam case, the initial intermediate images are processedhierarchically by the algorithm. Analogously to the fan-beam case,intermediate images that are close to each other in position andorientation are aggregated, in order that the aggregated intermediateimage might be sampled sparsely. The 3D sampling patterns for theintermediate images are determined by studying the structure ofconstituent weighted backprojected divergent-beams that are rotated andtranslated with respect to each other. One sampling pattern, asillustrated in FIG. 28, is where the samples are located along the raysof a divergent-beam corresponding to the central constituent beam of theintermediate image. The sample-spacing along each ray is chosen toensure that all the constituent divergent-beams of that intermediateimage are sufficiently sampled along the ray. Alternatively, a moregeneral way is to use the local Fourier method to find the samplingpattern, as described previously for the fan-beam case. Knowing how anintermediate image is composed of its constituent projections, the local3D Fourier structure of every intermediate image is determined. A 3Dlocal sampling matrix function at each point of the intermediate imageis found that matches the sampling requirements for the local Fouriersupport, as described in the fan-beam case. This matrix function is thenintegrated (possibly numerically) over the image domain to determine theposition of the samples.

A separable method of up-sampling combined with rotation and translationonto a new sampling beam is achieved similarly to the fan-beam case. Itreduces the 3D coordinate transformation and resampling operations to asequence of 1D resampling operations. As shown in FIG. 29(a), eachdivergent-beam may be regarded as the intersection of a set of verticalplanar fan-beams that are distributed in azimuthal angle, with a set oftilted planar fan-beams at different elevation angles. The steps of theseparable coordinate transformation are as follows (as shown in FIG.29(b))

-   -   1. The original divergent-beam is resampled onto the set of        intersection of the rays of the original with the vertical        planes of the new divergent-beam. These points are therefore        located on the planes shared by the final sampling points.    -   2. Steps 1 and 2 from the separable resampling in the fan-beam        case are performed for each plane separately to resample onto        the final set of points.

With a suitably efficient sampling scheme, the fast hierarchicalbackprojection algorithm for divergent-beam can be expected to achievelarge speedups with desirable accuracy. As in the fan-beam case onemight use a pseudo-beam that is modified (e.g., with the location of thevertex moving farther away from the origin) in successive levels.

As will be evident to those skilled in the art, the methods of thepresent invention are not limited to the examples of imaging geometriesor specific embodiments described herein. The methods are equallyapplicable to general problems of backprojection with other geometries.

FIG. 30 illustrates resampling-based backprojection, as the sum ofupsampled intermediate images formed by generalized backprojection ofindividual projections at source positions β_(p). It is similar to therotation-based backprojection 4, but differs from it in two respects.First, in the first step the p^(th) (possibly processed) projection issubjected to a weighted-backprojection 184 at the source position ororientation β_(p), rather than at zero, as is the case in 4. For examplein the case of fanbeam projections, each projection is backprojected atthe orientation of its source-angle. Second, before aggregation byaddition at 188, each initial intermediate image undergoes an upsamplingoperation at 186. This is necessary, because the sampling pattern ofeach of these P single-projection intermediate images is chosen to be anefficient and sparse pattern, so will usually be different for eachprojection, and often non-Cartesian. Before the intermediate images areaggregated, they need to be resampled onto a common and denser grid. Byitself, this structure offers no reduction in operation count. However,it serves as a stepping stone to the introduction of some of the fasthierarchical backprojection methods of the present invention.

FIG. 31 is another embodiment of the present invention, aresampling-based hierarchical backprojection for creating a pixel image{circumflex over (f)} from a plurality of projections q₁ . . . q_(p).FIG. 31 is the binary hierarchical version of FIG. 30. First, as in thenon-hierarchical case, each projection is backprojected 184 at itsindividual orientation, onto a sampling pattern suited to thatorientation, producing an intermediate digital image. This is the zerothlevel of the hierarchy. In the first level of the hierarchy, theseintermediate digital images are upsampled at 186 to a denser samplingpattern common to selected pairs of images. This upsampling will usuallybe to a non-Cartesian sampling pattern, but can be performed by asequence of one dimensional resampling operations, as shown in FIG. 23,or in FIG. 29(b). Selected resulting upsampled images are aggregatedpairwise at 190, producing new intermediate images. In the second level,the new intermediate images are again upsampled at 192, and aggregatedat 194, producing new intermediate images. This process continues untilall intermediate image and projections have been processed, producingafter the last aggregation step 198 the final image {circumflex over(f)}. As in the previously described embodiments, operations can becombined within and across a level.

The sampling patterns at each stage in the hierarchy and the parametersof the upsampling operation in the embodiment shown in FIG. 31 may bechosen by any of the previously described methods. For example, in thecase of fan-beam projections, one possible sampling pattern for a givenintermediate image would lie on the points of intersections of two fans:the first oriented at the central constituent source position; thesecond oriented at an extremal constituent source position.Alternatively the sampling pattern and parameters of the upsamplingsteps can be determined based on the Fourier or local Fourier propertiesand the view angle of projections included in the intermediate images.In particular, for non-periodic sampling patterns, the local Fouriermethod described for the fanbeam rotation-based algorithm can be used tofind the sampling patterns: knowledge of how the projections combine toform the intermediate images leads to the determination of the localspectral support, which is used in turn to calculate the local samplingmatrix function, which when integrated over the image domain producesthe sampling pattern for the intermediate image.

FIG. 32 is the block-diagram for fast hierarchical reprojection.Reprojection is the process of computing tomographic projections from agiven electronic image. The reprojection algorithm is found by applyingthe process of flow graph transposition to any block diagram of a abackprojection algorithm, possibly with some change in weighingoperations. In the process, operations are replaced by their adjoint ordual. The block diagrams for reprojection therefore appear similar to areversed version of the corresponding one for backprojection. Thereprojection process described in FIG. 32 is one such embodiment ofreprojection, obtained from the backprojection algorithm described inFIG. 8.

In the first level, a copy of the input image ƒ is preserved in the topbranch of the diagram as a top intermediate image, and in the bottombranch the image ƒ is rotated at 200 by π/2, producing a bottomintermediate image. In the second level, in the top half of the diagram,the un-rotated to intermediate image is subject to three separatedigital image coordinate transformations at 202 some of which may leavethe image unchanged, producing three different top intermediate images.A similar process is applied in the bottom branch, producing threebottom intermediate images. Each of the top and bottom intermediateimages (six in all) then undergoes a process of decimation (low-passfiltering followed by downsampling) at 204, producing new intermediateimages. In the instance of the embodiment illustrated in FIG. 32 thereare only 2*3^(L), with L=2 view-angles, so the third level is the finalone in the recursive hierarchy. In the third and final level theintermediate images are subject to separate coordinate transformations(206) some of which may leave the image unchanged, producing 18intermediate images. The last step, which is not part of the recursionis different: each intermediate image undergoes a reprojection at zerodegree at 208. Reprojection at zero degrees is equivalent to summing thevertical lines of pixels to produce a one-dimensional signal. Theseone-dimensional signals (518) are the output projections produced by thealgorithm.

The parameters of the digital image coordinate transformations in thealgorithm are chosen by the knowledge of the Fourier characteristics ofthe intermediate images. These parameters are simply related to theparameters of the corresponding backprojection algorithm. It is easy tosee that since the reprojection block-diagram is a flow transposition ofa backprojection block-diagram, every branch of the reprojectionblock-diagram has a corresponding branch in the backprojectionblock-diagram, and the coordinate transformations in the correspondingbranches are mathematical adjoints of each other. In the version of thisreprojection algorithm that corresponds to the two-shear backprojectionalgorithm in FIG. 10, the coordinate transformations in the second level(202) is an x-shear followed by a y-shear, and the coordinatetransformation in the last level (206) is an x-shear,followed by ay-shear, and a fractional decimation in x (These three operations can bereduced to a shear-scale). The parameters of these shears are thenegative of the corresponding parameter used in the backprojectionalgorithm. The parameter of the decimation in x is the same as that ofthe upsampling in x in the first level of the backprojection algorithm.In the shear-scale version of this algorithm (corresponding to theshear-scale backprojection algorithm displayed in FIG. 13), thecoordinate transformations in the second level (202) are shears in x(and the parameter of each shear is the negative of the correspondingparameter in FIG. 13). The coordinate transformations in the final level(206) are shear-scales. The shear-scale used in the reprojectionalgorithm is the mathematical adjoint of the corresponding shear-scaleused in the backprojection. The parameters of the decimation factors arealso the same as the corresponding upsampling factors in thebackprojection.

Just like in the backprojection algorithms, oversampling of theintermediate images in the algorithm can be enforced by first upsamplingthe images at the beginning of the algorithm and downsampling them bythe same factor at the end of the algorithm. Also, these operations canbe combined within or across a level.

FIG. 33 is the block-diagram of a decimation-based weightedreprojection. It it the flow-graph transpose of FIG. 31 It shows thereprojection of the image onto a set of projections at 18 differentsource angles. Initially the given image is processed along two parallelpaths. In each path the image is subject to three parallel resamplings(210) onto three different sparser sampling patterns. This resamplingcan be performed in a separable way using one-dimensional decimations(low-pass filtering followed by resampling onto a sparser grid). Theparameters of the filter are determined by the Fourier characteristicsof the intermediate image and the desired projections. Local Fourieranalysis of the desired projections is used in the case when projectionsdo not line on parallel straight lines. In the final level of thealgorithm each intermediate image is again subjected to three parallelresamplings (212) onto sparser sampling patterns. Finally a weightedprojection is performed on the image to produce the projection p_(θ).This involves a weighted sum of the pixels of the image to produce aone-dimensional projection.

Implementations and Experimental Results

Preferred embodiments of the present invention were implemented in Cprogramming language and tested in numerical experiments on a Sun Ultra5 workstation with 384 MB RAM. The test image was the Shepp-Logan headphantom(a standard test image used in numerical evaluations oftomographic algorithms) By varying the parameters of the algorithms atradeoff can be made between accuracy and computational cost. Accuracyrefers to the quality of the reconstructed image. Though visual qualityis not easily quantifiable, we measure the error between thereconstructed image and the original from which the radon transform wasnumerically computed. The measure of error used is the relativeroot-mean-square error (RRMSE). The RRMSE in reconstructing an N×N image{circumflex over (f)}[m₂, m₁] from the tomographic projections of f[m₂,m₁] is calculated as follows: $\begin{matrix}{{RRMSE} = \frac{\sqrt{\frac{1}{N^{2}}{\sum\limits_{m_{2} = 1}^{N}{\sum\limits_{m_{1} = 1}^{N}{{{\hat{f}\left\lbrack {m_{1},m_{2}} \right\rbrack} - {f\left\lbrack {m_{1},m_{2}} \right\rbrack}}}^{2}}}}}{\max_{m_{1},m_{2}}{f\left\lbrack {m_{1},m_{2}} \right\rbrack}}} & (1)\end{matrix}$

For parallel-beam data, the test image was of size 256×256, the numberof view angles was 486, and the Shepp-Logan filter (the ideal rampfilter multiplied by a sinc function) was used to filter individualprojections. FIG. 30 displays the RRMSE error versus the run times forthe two algorithms at various values of the oversampling parameter γbetween 0.75 and 1.0. The two-shear algorithm is represented by thecircles and the shear-scale algorithm by the squares. Each algorithm isrun using two types of filters—a third-order (dashed line) andfifth-order (solid line) spline filter called MOMS 16. For each flavorof the algorithm, as γ is decreased the error of the algorithm decreasesand the run time increases. The plot points that are not connected toany other are the non-oversampled versions of the algorithms representedby the connected points. In comparison the run time of the conventionalalgorithm, using linear interpolation, is 14s and the RRMSE error of itsoutput image is 0.0486 (worse than the fast algorithms displayed here).

Some sample images from the output of the algorithms, for parallel-beamdata, are displayed in FIG. 35. Columnwise from left to right, they areoutput images from the conventional backprojection, the two-shear andthe shear-scale algorithms. An oversampling of γ=0.82 was applied to thetwo fast algorithms. The lower row of images displays in detail asection of the corresponding images in the upper row.

The fast hierarchical algorithm for fan-beam geometry was successfullytested on the 512×512 2D shepp-logan phantom. The acquisition geometryconsidered was with a source radius D=1.06×N=544, 972 source angles, and1025 equiangular detectors. The regular and oversampled version of thefast algorithm was implemented using a variety of interpolation methods.The resulting reconstructed images were compared to a conventionalfan-beam algorithm that used linear interpolation. In all theexperiments the projections were filtered with the Shepp-Logan filter.

FIGS. 36(A) and 36(B) display the reconstructed images from theconventional in FIG. 36(A) and the fast algorithms in FIG. 36(B). Theresult of the fast algorithm is comparable to that of the conventionalalgorithm. The point spread functions of the fast algorithms arecomparable to that of the conventional one. FIG. 37(B) displays the PSFof the conventional algorithm, FIG. 37(C) displays the PSF of the fastalgorithm with linear interpolation and γ=0.4 oversampling. FIG. 37(c)compares slices through the x-axis of the psfs of the conventionalalgorithm, the non-oversampled and the γ=0.4 oversampled fastalgorithms. The similarity of the PSFs confirms the comparable imagequality in the fast algorithms.

The sampling scheme used was the fan intersection-based method, withoutthe enhancements suggested later. In addition to the shortcomings ofthis scheme mentioned previously, for reasons of simplicity ofimplementation, numerous sample points not needed for the correctoperation of the algorithm were used in the embodiment tested in theexperiments. Despite these inefficiencies, for N=512 and D=1.41N, theratio of (data-dependant) multiply operations in the conventionalalgorithm to the fast (linear) one was 6.4 and the ratio of additionoperations was 3.0. Note that the geometry used in this experiment, witha source very close to the origin (D=2.06R) is particularly challengingfor this un-enhanced implementation, because of the high density ofsamples near the source vertex. In most practical systems, the source isfurther away, reducing these effects. Furthermore using the alternativesampling schemes discussed earlier, will lead to much higher speedupfactors.

Detailed Description of Backprojection and Reprojection Algorithms Usedin the Present Invention

Backprojection is the process used to create images from processedprojections. Reprojection is the reverse process, used to computeprojections from a given image. Both operations are used in imagereconstruction from projections, as seen in FIG. 1. Conventionalbackprojection will be described first, followed by a description of thebackprojection and reprojection methods of the present invention. Thedescription will be given first for the case of parallel-beamprojections of a 2D image, and then for the more general 2D and 3Dgeometries.

The classic and preferred method used to estimate an image from itsprojections is the filtered backprojection (FBP) algorithm. The FBPinvolves first filtering the projections with a so-called ramp ormodified ramp filter and then backprojecting those filtered projections.Additional pre-processing steps may also be applied to the projectionsbefore backprojection. A processed projection taken at view angle θ_(p),where p=1, . . . P is the projection index, will be denoted alternatelyby q(t, p), or q_(p) or q_(θ) _(p) , depending on whether the dependenceon the variable t or θ needs to be shown explicitly. For brevity,processed projections will be called projections hereinafter.

The FBP reconstruction {circumflex over (f)} from a set of P projections(at the angles specified by the components of the length-P vector {rightarrow over (θ)}), is described by:{circumflex over (f)}(x ₁ ,x ₂)=(B _({right arrow over (θ)}) q)(x ₁ ,x₂)   (2)where

Definition 0.1 the backprojection operator is defined by $\begin{matrix}{{\left( {\mathcal{B}_{\overset{\rightarrow}{\theta}}q} \right)\left( {x_{1},x_{2}} \right)}\overset{\bigtriangleup}{=}{\sum\limits_{p = 1}^{P}{q\left( {{{x_{1}\cos\quad\theta_{p}} + {x_{2}\sin\quad\theta_{p}}},p} \right)}}} & (3)\end{matrix}$

The O(N² log P) coordinate-transformation-based hierarchicalbackprojection algorithms of the present invention are based on ahierarchical decomposition of backprojection in terms of coordinatetransformations. The details of several preferred embodiments aredescribed with different choices of coordinate transformation,concluding with the most general coordinate transformations.

In particular, the O(N² log P) rotation-based hierarchicalbackprojection algorithms of the present invention are based ondecomposition of backprojection in terms of the rotation of images.

The rotation-by-θ operator K(θ) which maps an image ƒ( {right arrow over(x)}) to its rotated version (K(θ)f)({right arrow over (x)}) is definedby

Definition 0.2 (K(θ)f)({right arrow over (x)})

ƒ( K_(θ){right arrow over (x)}),

where

Definition 0.3 The matrix of rotation by angle θ in the plane is$K_{\theta} = \begin{bmatrix}{\cos\quad\theta} & {{- \sin}\quad\theta} \\{\sin\quad\theta} & {\cos\quad\theta}\end{bmatrix}$

Backprojection of a 1D function q(t) at a single-angle θ is defined as(B_(θ)q)({right arrow over (x)})

πq(x₁ cos θ+x₂ sin θ). In particular, the zero-angle backprojection (forθ=0) is(B ₀ q)({right arrow over (x)})=q(x ₁ cos 0+x ₂ sin 0)=q(x ₁)=q([10]{right arrow over (x)})   (4)

Definition 0.4 The single-θ_(p)-backprojection intermediate image isI_(θ) _(p) ({right arrow over (x)})

(B_(θ) _(p) q_(p))({right arrow over (x)})

As seen in FIG. 4, the backprojection in equation (3) may be rewrittenas follows. Here q_(p) refers to the projection whose view-angle isθ_(p). $\begin{matrix}{{\hat{f}\left( \overset{\rightarrow}{x} \right)} = {\sum\limits_{p = 1}^{P}{\left( {{\mathcal{K}\left( {- \theta_{p}} \right)}\mathcal{B}_{0}q_{p}} \right)\left( \overset{\rightarrow}{x} \right)}}} & (5)\end{matrix}$

One embodiment of hierarchical backprojection stems from the fact thatthe cumulative result of several successive rotations is still arotation. In particular,K(θ₁)K(θ₂) . . . K(θ_(N))=K(θ₁+θ₂+ . . . +θ_(N))   (6)

It follows from Equation (6) that with P=2^(L) for some integer L, theblock diagram in FIG. 4 can be rearranged into a hierarchical treestructure as shown in FIG. 5. The intermediate image in the m^(th)branch of the l^(th) level is denoted as I_(l,m). In the initial levelI_(l,m)

B₀q_(m) m=1, 2, . . . , P   (7)intermediate images are produced by backprojecting individualprojections, and in subsequent levels, i.e., for I=1, 2, 3, . . . , log₂PI_(l+1,m)

K(δ_(1,2m−1))I_(l,2m−1)+K(δ_(l,2m))I_(l,2m)   (8)the intermediate images undergo coordinate transformations (rotations,in this embodiment) and are aggregated by addition. For any set ofprojection angles θ_(i), i=1, . . . , P, the intermediate rotationangles δ_(l,m), can be chosen to guarantee that the structures in FIG. 4and 5 are equivalent, so that the hierarchical backprojection algorithmdepicted in FIG. 5 produces the desired image {circumflex over (f)}.

In fact, because there are many more intermediate rotation angles (freeparameters) than view angles (constraints), there are many degrees offreedom in the choice of the δ_(l,m). Also, for the digital images usedin practice, digital image coordinate transformations are used, whoseaccuracy and computational cost depend on the choice of samplingpatterns and implementation. These various degrees of freedom, orparameters, are used in the present invention to reduce thecomputational requirements of the hierarchical backprojection algorithm,as will be described.

Definition 0.5 The set N_(l,m) of indices of constituent projections foran intermediate image I_(l,m) is the set of indices of projections forwhich there is a path from the input in FIG. 5 to the intermediate imageI_(l,m).

For example, as seen in FIG. 5, N_(3,P/4)={P−3, P−2, P−1, P}. It iseasily shown, using equation (8) or upon examination of the blockdiagram, thatN _(l,m)={2^(l−1)(m−1)+k: k=1, 2, . . . , 2^(l−1)}  (9)

Thus N_(l,m) lists the projections that make up intermediate imageI_(l,m). If, instead of being processed by the method described in FIG.5, these same projections indexed by the set N_(l,m) were directlybackprojected together at their respective view-angles, this wouldproduce the virtual intermediate image Ĩ_(l,m):

Definition 0.6 Ĩ_(l,m)=Σ_(pεN) _(l,m) I_(θ) _(p)

Upon examination of the block diagram in FIG. 5, it can be seen that therelative angle between projections in an intermediate image is preservedin the final reconstructed image {circumflex over (f)}. In other words,for all intermediate images in the algorithm, i.e. ∀l, mI _(l,m) =K(α_(l,m))Ĩ _(l,m) for some α_(l,m) ε(−π, π]  (10)

The intermediate rotation angles δ_(l,m) of the hierarchical algorithmare completely determined by α_(l,m) as follows. It follows easily fromthe definition of N_(l,m) or Equation (9) that N_(l+1,m)=N_(1,2m−1) ∪N_(l,2m), and consequently by the Definition Ĩ.Ĩ _(l+1,m) =Ĩ _(l,2m−1) +Ĩ _(1,2m)   (11)Equations (8), (10) and (11) imply thatδ_(l,2m−1)=α_(l+1,m)−α_(l,2m−1) and δ_(l,2m)=α_(l+1,m)−α_(l,2m)   (12)

The intermediate rotation angles (δ_(l,m)) and correspondingly therelative angles α_(l,m), are chosen to reduce the bandwidth of theintermediate images. Images of small bandwidth can be represented by fewsamples, and consequently reduce the computational cost of the wholealgorithm. The bandwidth of these intermediate images is determined byunderstanding the algorithm in the Fourier domain and tracking theevolution of the spectral support of the intermediate images through thehierarchy.

Tomographic Projections in the Fourier Domain

The parameters of the digital image coordinate transformations arechosen to account for the view angles θ_(p) of the selected projections,as has been described. The parameters are also chosen for their effecton the Fourier characteristics of the intermediate images. Theseconsiderations are also used to determine which intermediate images areselected for aggregation together.

Key to interpreting the backprojection algorithm in the Fourier domainis the projection-slice theorem that relates the one-dimensional Fouriertransform (F₁) of a projection P_(θ)f with the two-dimensional Fouriertransform (F₂) of the image ƒ. The projection-slice theorem 4 says that(F ₁ P _(θ) f)(ω)=(F ₂ f)(ω cos θ, ω sin θ)

The backprojection algorithm of FIG. 4 can be interpreted in the Fourierdomain. The backprojection-at-zero operator produces an image whosespectral support is limited to the horizontal frequency axis ω₁:(B ₀ q)(x ₁ ,x ₂)=q(x ₁)

(F ₂ B ₀ q)(ω₁,ω₂)=2πQ(ω₁)δ(ω₂)   (13)

It is well known that rotating an image in the space domain results inthe rotation of its Fourier transform by the same angleg({right arrow over (x)})=(B ₀ q)(K _(θ) {right arrow over (x)})

G({right arrow over (ω)})=(F ₂ B ₀ q)(K _(θ{right arrow over (w)}))  (14)

It follows that the function of the backprojection operator in Equation(5) is to rotate the spectral component of each projection to theappropriate angle in the reconstructed image and add them all together.

Assuming projections with one-sided bandwidth (in the t variable) equalto π, the typical virtual continuous intermediate image Ĩ_(l,m) in thehierarchical algorithm therefore has a spectral support of a wedge in arange of angles determined by the constituent projections, as shown inFIG. 6. In particular, let I_(l,m) be an intermediate image in the blockdiagram of the algorithm shown in FIG. 5(A), with N_(l,m)={b, b +1, . .. , e}, i.e. the view-angles of the constituent projections of thisimage are {θ_(p): p=b, b +1, . . . , e}. Because of the relationship inEquation (10), the spectral support of I_(l,m) is just a rotated versionof that of Ĩ_(l,m), with rotation angle α_(l,m). It follows that therange of angles occupied by the wedge in FIG. 6 is [γ−β,γ+β]=[θ_(b)−α_(l,m), θ_(e)−α_(l,m)]. The bandwidths of I_(l,m) in thefirst and second coordinate, Ω_(s) and Ω_(f), respectively, areindicated on the bounding rectangle for the wedge, which is shown bybroken lines in FIG. 6. Thus, the choice of α_(l,m) provides a way tocontrol the bandwidths of I_(l,m).

Sampling theory dictates how such a continuous intermediate imageI_(l,m) may be represented by a samples on a Cartesian sampling pattern.In particular, a continuous image with bandwidths of Ω_(s) and Ω_(ƒ) inthe first and second coordinate may be sampled on a Cartesian patternaligned with the coordinate axes, with a sample spacing of Δ_(s) andΔ_(ƒ) in the first and second coordinates such that Δ_(s)<π/Ω_(s) andΔ_(f)<π/Ω_(f). In order to minimize the number of samples required torepresent the continuous image, the quantity1/(Δ_(s)Δ_(f))>π²Ω_(s)Ω_(f), which is proportional to the area of thebounding rectangle, should be minimized. The rotation angle thatminimizes this area, and therefore minimizes the sampling requirementsfor intermediate image I_(l,m), isα_(l,m)*=(θ_(b)+θ_(e))/2   (15)Selecting Parameters for Coordinate Transformations

Equation (15) is used in Equation (12) to determine the optimum rotationparameters δ_(l,m) for the digital image coordinate rotations in therotation-based hierarchical backprojection embodiment illustrated inFIG. 5. The indices b and e are determined as corresponding to thesmallest and largest, respectively, in the set N_(l,m). With thischoice, the intermediate image I_(l,m)=K(α_(l,m)*)Ĩ_(l,m), and itsbandwidths in the first coordinate (slow bandwidth) and secondcoordinate fast bandwidth) are, respectively, $\begin{matrix}{{{\Omega_{s}\left( I_{l,m} \right)} = {\pi\quad\sin\quad\frac{\theta_{e} - \theta_{b}}{2}}}{{\Omega_{f}\left( I_{l,m} \right)} = \pi}} & (16)\end{matrix}$

To further minimize the bandwidth, the differences θ_(e)θ_(b) should beminimized, which is achieved by selecting to aggregate at each stepintermediate images produced from projections with the least maximumangular separation between their view angles.

In addition to rotation angles and selection of which intermediateimages to aggregate, it is necessary to choose appropriate samplingpatterns for the digital image coordinate transformations of the variousintermediate images in the algorithm. These are chosen using thesampling requirements, which are in turn determined using the spectralsupports and bandwidths of the intermediate images.

In particular, for initial intermediate images formed by thebackprojection of a single projection, I_(l,m)=B₀q_(θ) _(m) as in FIG.5, the slow bandwidth Ω_(s)(I_(l,m))=0, and a single sample along the x₂coordinate suffices. The slow bandwidth will be larger and more samplesalong the x₂ coordinate will be needed for initial intermediate imagesformed, as in FIG. 3-b, by backprojecting more than one projection.

The sum of two virtual intermediate images(Ĩ_(l,m)=Ĩ_(l−1,2m−1)+Ĩ_(l−1,2m)) and their corresponding intermediateimages are illustrated, in the space and frequency domains, in FIGS.7(A)-7(C). FIGS. 7(A) and 7(B) show the Fourier-domain support of Ĩ andI. TFig. 7(C) shows the space-domain sampling scheme ofI_(l,m),I_(−1,2m) and I_(1−1,2m−1). The sampling points are at theintersections of the horizontal and vertical lines in the space domain.The aggregated intermediate image I_(l,m) having a greater slowbandwidth than each of its components, requires a denser samplingpattern. Conversely, intermediate images at earlier levels of thehierarchical algorithm require a sparser sampling pattern, leading tocomputational savings.

Computation Cost and Savings

The computational cost is readily estimated. Assuming equally Pprojection view-angles in [0, π), with an angular spacing of Δ_(θ)=π/P,and N_(l,m) as specified in (9), equation (16) (using b=2⁻¹(m−1)+1 ande=2^(l−1)(m−1)+2^(l−1)) implies that Ω_(s)(I_(l,m))=π sin(π(2⁻¹−1)/(2P))and Ω_(f)(I_(l,m))=π. The size (number of points) of the digital imageI_(l,m) ^(d) in the l^(th) level is therefore $\begin{matrix}{{{size}\left( I_{l,m}^{d} \right)} \approx {\left( {N/\Delta_{s}} \right)\left( {N/\Delta_{f}} \right)}} \\{= {\left( {N\quad{{\Omega_{s}\left( I_{l,m} \right)}/\pi}} \right)\left( {N\quad{{\Omega_{f}\left( I_{l,m} \right)}/\pi}} \right)}} \\{= {{N^{2}{\sin\left( {{\pi\left( {2^{l - 1} - 1} \right)}/\left( {2P} \right)} \right)}} < {N^{2}\left( {{\pi 2}^{l - 1}/\left( {2P} \right)} \right)}}} \\{= {{O\left( {N^{2}{2^{l}/P}} \right)}.}}\end{matrix}$

Given that the cost of rotating a digital image of S pixels is O(S) andsize(I_(l,m) ^(d))=O(N²2^(l)/P), the complexity of the algorithm is$\begin{matrix}{{\sum\limits_{l = 1}^{\log\quad P}{\frac{P}{2^{l - 1}}{O\left( \frac{N^{2}2^{l}}{P} \right)}}} = {{\sum\limits_{l = 1}^{\log\quad P}{O\left( N^{2} \right)}} = {O\left( {N^{2}\log\quad P} \right)}}} & (17)\end{matrix}$

For the typical P=O(N), this is O(N² log N), which is much morefavorable scaling of the cost with the size of the image, than the O(N³)of conventional backprojection.

Improved Ternary Rotation-Based Hierarchical Backprojection

It is easy to see that the rotation of a sampled image by certainspecial angles—0,±π/2 and π—are computationally inexpensive operationsbecause they involve no interpolation but, at most, a mere rearrangingof existing pixels. The computational efficiency of the algorithm of thepresent invention can be improved by incorporating these free operationsinto it. The binary algorithm (FIG. 5) can be modified such that three,not two, images are added at every stage. The center image of each suchtriplet is rotated by 0 radians—a free operation. To use the freerotation by −π/2, it is included in the final stage: one of theconstituent images is rotated by −π/2 and added to the other un-rotatedone. This particular combination of binary and ternary stages results ina set of view-angles of size 2×3^(L) for some integer L. Though thealgorithm can be tailored to arbitrary sets and numbers of projections,the description and analysis of the algorithm can be simplified byassuming that they number exactly 2×3^(L) and are uniformly distributedas follows: $\begin{matrix}{{\theta_{i} = {{{- \frac{\pi}{4}}\left( {1 - 3^{- L}} \right)} + {{\Delta\quad}_{\theta}\left( {i - 1} \right)}}}{where}{{i = 1},2,\ldots\quad,{{{2 \cdot 3^{L}}{and}\Delta_{\theta}} = {\pi/\left( {2 \cdot 3^{L}} \right)}}}} & (18)\end{matrix}$

Hence, the view-angles can be divided into two sets of 3^(L) each, oneset centered around θ=0 and the other centered around θ=π/2. The blockdiagram of this ternary algorithm is shown in FIG. 8 for the particularcase of L=2, i.e., for a set of P=2×3²=18 projections. The digitalintermediate images I_(l,m) ^(d) are the sampled versions of theunderlying continuous intermediate images I_(l,m). In the block diagram,blocks marked B₀ represent the zero-angle backprojection operator.Blocks marked K_(δ) represent the digital image rotation operator. Aswill be explained, the sampling periods and rotation angles for thedigital images and for the digital image rotations are chosen tominimize computational cost.

Selected aggregate intermediate digital images can be stretched toproduce new intermediate images. The blocks labelled ↑_(y) U_(l,m) inFIG. 8 represent up-sampling along the second (slow) coordinate, whichare necessary to accommodate the increasing bandwidth of theintermediate images as the algorithm progresses. For the configurationof view-angles in (18) and (22), the bandwidth of the intermediateimages in the l^(th) level are Ω_(s)(I_(l,m))=π sin(Δ_(θ)(3^(l−1)−1)/2)and Ω_(f)(I_(l,m))=π. The normalized spectral supports at key points ofthe algorithm, numbered (1), . . . , (9), are displayed in FIG. 9. Asthe algorithm progresses (as l increases), the slow bandwidth increases,and the required slow-sampling period,Δ_(s)<π/Ω_(s)(I_(l,m))=sin((3^(l−1)−1)/2)⁻¹, decreases. The upsamplingfactors U_(l,m) in the various upsampling blocks can be chosen tosatisfy these requirements.

Section 0.0.0.1 describes how the separable rotation of discrete imagesinvolving two-shears is incorporated into the backprojection algorithm.In order to improve the quality of the backprojected image, a systematicoversampling of the intermediate images is introduced. Section 0.0.0.1describes the algorithms modified to include this oversampling.

Intermediate rotation angles for coordinate transformation can be chosenas follows. For general L, the hierarchy consists of (L+1) levels. Inthe zeroth level, each projection q_(θ) is backprojected-at-zero (B₀) toproduce a corresponding image B₀q_(θ). The images are then grouped intothrees and combined to produce a third as many images. The groupings atlevel l are defined by the following relationships. In subsequentlevels, l=1, . . . , L,I _(l+1,m) =K(δ_(l,3m−2))I _(l,3m−2) +K(δ_(l,3m−1))I _(l,3m−1)+K(δ_(l,3m))I _(l,3m)   (19)and I _(L+2,1) =K(δ_(L+1,1))I _(L+1,1) +K(δ _(L+1,2))I _(L+1,2)   (20)

For the configuration of view angles in (18), the optimal intermediaterotation angles are δ_(l,3m−1)=0 for l=1, 2, . . . , L , δ_(L+1,1)=0,and δ_(L+1,2)=−π/2. Upon examination of FIG. 10, it is easy to see thatthe indices N_(l,m) of the projections that constitute image I_(l,m) areas follows:N _(l,m)={3^(l−1)(m−1)+1, 3^(l−1)(m−1)+2, . . . , 3^((l−1)) m} for l=1,2, . . . , L+1   (21)

By Equation (21) therefore, b=3^(l−1)(m−1)+1 and e=3^(l−1)(m−1)+3^(l−1)and, using Equations(18) and (15)α_(l,m)*=−π/4(1−3^(−L))+Δ_(θ)(3^(l−1)(m−½)−½). Similarly to the binaryit follows that for l=1, . . . L, $\begin{matrix}{\delta_{l,{{3m} - i}} = {{\alpha_{{l + 1},m}^{*} - \alpha_{l,{{3m} - i}}^{*}} = \left\{ \begin{matrix}{{+ \Delta_{\theta}}3^{l - 1}} & {{{if}\quad i} = 2} \\0 & {{{if}\quad i} = 1} \\{{- \Delta_{\theta}}3^{l - 1}} & {{{if}\quad i} = 0}\end{matrix} \right.}} & (22)\end{matrix}$

For the final stage we make use of the free rotations by π/2 and 0radians.

0.0.0.1 Coordinate transformation can be based on shearing, as follows.The digital image rotation of the intermediate images can be replaced bya sequence of two digital image shears as shown in FIG. 10(B), whereshears in the x and y coordinate are defined as follows

Definition 0.7 The x-shear and y-shear operators are defined by(S _(x)(α)f)({right arrow over (x)})=f(S _(x) ^(α) {right arrow over(x)})(S _(y)(α)f)({right arrow over (x)})=f(S _(y) ^(α) {right arrow over(x)})where the x-shear and y-shear matrices are${S_{z}^{\alpha} = {{\begin{bmatrix}1 & \alpha \\0 & 1\end{bmatrix}\quad{and}\quad S_{y}^{\alpha}} = \begin{bmatrix}1 & 0 \\\alpha & 1\end{bmatrix}}},$respectively.

This alternative embodiment is derived from the two-shear factorizationof the rotation matrix: K_(θ)=S_(y) ^(tan θ)S_(x) ^(−sin θcos θ)S_(c),where $S_{c} = {\begin{bmatrix}{\cos\quad\theta} & 0 \\0 & {{1/\cos}\quad\theta}\end{bmatrix}.}$This decomposition implies that in the case of digital images, usingperfect bandlimited interpolation and a sufficiently bandlimited digitalimage as the input to the two-shear digital image coordinatetransformation, the twice-sheared image is simply the image rotated byθ, and effectively downsampled in x and upsampled in y by a factor of1/cos θ. To correct these changes in sampling pattern would requirenon-integer resampling every time a two-shear operation is performed.Instead, by sampling the intermediate images differently, only thecumulative effects need to be addressed for the final image.

A convenient choice is to correct the resampling once at the beginningof the hierarchy, by upsampling the initial digital intermediate imagesin x and downsampling in y. When the initial intermediate images areproduced by single-view backprojection as in the embodiments shown inFIG. 10, the downsampling in y is free, because it only involves initialbackprojection onto a sampling pattern with a different y density. It istherefore not shown explicitly in the block diagram in FIG. 10.Upsampling in x in the first stage avoids aliasing problems in laterstages. For each of the initial images, the cumulative effectivedownsampling/upsampling factor 1/Π_(k=1) ^(K) cos θ_(k) to the finalimage is calculated across all two-shear transformations on the path tothe final image, and the reverse operation is applied at the initialimage as described. It follows that each of the initial images I_(1,p)^(d), p=1, . . . , P is resampled with different x and y samplingdensities.

So the resulting algorithm, as shown in FIG. 10(A) involves firstresampling all the images B₀q_(θ) _(p) , p=1, . . . , P by theappropriate factors (depending on θ_(p)) and then performing two-sheardigital image coordinate transformations as shown in FIG. 10(B). The twoshears are performed taking into consideration the changing samplingpattern of each intermediate image.

Shear-Scale Algorithms

The shear-scale-based hierarchical backprojection algorithm is based onthe definition of backprojection as the sum of single-projection imagesthat are scaled (stretched/contracted) and sheared.

Definition 0.8 The x-shear-scale operator C({right arrow over (σ)}) isdefined by(C({right arrow over (σ)})f)({right arrow over (x)})=f(C_({right arrow over (σ)}) {right arrow over (x)})where the x-shear-scale matrix$C_{\overset{\rightarrow}{\sigma}} = {\begin{bmatrix}\sigma_{1} & \sigma_{2} \\0 & 1\end{bmatrix}.}$

So, Equation (3) may be rewritten as $\begin{matrix}{{\hat{f}\left( \overset{\rightarrow}{x} \right)} = {\sum\limits_{p = 1}^{P}{\left\lbrack {{\mathcal{C}\left( {{\cos\quad\theta_{p}},{\sin\quad\theta_{p}}} \right)}\mathcal{B}_{0}q_{p}} \right\rbrack\left( \overset{\rightarrow}{x} \right)}}} & (23)\end{matrix}$

The following are some results describing the effect of accumulatedshears and shear-scales that can be easily shown. $\begin{matrix}{{\prod\limits_{i = 1}^{n}\quad{S_{x}\left( \alpha_{i} \right)}} = {S_{x}\left( {\sum\limits_{i = 1}^{n}\alpha_{i}} \right)}} & (24) \\{{{\mathcal{C}\left( {\overset{\rightarrow}{\sigma}}^{\prime} \right)}{\mathcal{C}\left( {\overset{\rightarrow}{\sigma}}^{\prime\prime} \right)}} = {\mathcal{C}\left( {{\sigma_{1}^{\prime}\sigma_{1}^{\prime\prime}},{\sigma_{1}^{\prime\prime} + {\sigma_{1}^{\prime}\sigma_{2}^{\prime\prime}}}} \right)}} & (25)\end{matrix}$

Similarly to Equation (24), Equation (25) can be inductively applied toprove that a sequence of n shear-scales is also a shear-scale. Usingthis property, Equation (23) can be rewritten in a hierarchicalalgorithmic structure. FIG. 11(A) is the block-diagram representation ofEquation(23) (for P=4), and FIG. 11(B) is its hierarchical equivalent.

The requirement of equivalence of the two structures sets up a system ofequations between the given projection-angles θ_(p) and the unknownshear-scale parameters {right arrow over (σ)}_(l,m), whose solutionalways exists.

It is computationally beneficial (in terms of operations or storagespace) to use an integer scale-factor rather than a non-integer one. Forexample, it is beneficial to use x-shear-scales C({right arrow over(σ)}) with an integer shear-parameter σ₁=1.0, which is exactly a purex-shear S_(x)(σ₂). It is possible to design the algorithm so that onlypure x-shears are used, and x-shear-scalings are avoided, in levels l=2,3, . . . , L.

Instead of replacing each shear-scale in Equation (23) with a sequenceof shear-scales, using Equations (24) and (25) each shear-scale isreplaced with a single shear-scale followed by a sequence of shears:$\begin{matrix}\begin{matrix}{{\mathcal{C}\left( {{\cos\quad\theta_{p}},{\sin\quad\theta_{p}}} \right)} = {{\mathcal{S}_{x}\left( \alpha_{n} \right)}{\mathcal{S}_{x}\left( \alpha_{n - 1} \right)}\ldots\quad{\mathcal{S}_{x}\left( \alpha_{1} \right)}{\mathcal{C}\left( {{\cos\quad\theta_{p}},\sigma} \right)}}} \\{{{s.t.\quad\sigma} + {\cos\quad\theta_{p}{\sum\limits_{i = 1}^{n}\quad\alpha_{i}}}} = {\sin\quad\theta_{p}}}\end{matrix} & (26)\end{matrix}$with free parameters α₁, α₂, . . . , α_(n).

The resulting hierarchical configuration is displayed in FIG. 13. Aconvenient choice is to perform the x-shear-scale(C_({right arrow over (σ)})) in the level l=1 and x-shears in subsequentlevels l=2, 3, . . . , L. So the expression for the underlyingcontinuous-domain image in the first level is as follows.$\begin{matrix}{I_{2,m} = {\sum\limits_{i = 0}^{2}\quad{\mathcal{C}{\overset{->}{\left( \sigma^{{3m} - i} \right)\mathcal{B}_{0}q_{{3m} - i}}}^{\quad}}}} & (27)\end{matrix}$For l=2, 3, . . . , L the intermediate image I_(l+1,m) is composed ofthe x-sheared versions of three intermediate images ({I_(l,3m−1)}_(i=0)²) in the previous level:I _(l+1,m) =S _(x)(α_(l,3m−2))I _(l,3m−2) +I _(l,3m−1) +S_(x)(α_(l,3m))I _(l,3m)   (28)And in order to use the free rotation by 0 and −π/2 radians, in thefinal level: I_(L+2,1)=I_(L+1,1)+K(−π/²)I_(L+1,2).

The parameters of the coordinate transformations in the algorithmdisplayed in FIG. 13 (for the various shear-scales and shears) areunderdetermined. In fact, only the (nonintegral) scale-factors σ₁ in thefirst level are completely determined by the angle of the associatedprojection. The intermediate shear-factors({α_(l,m)}_(l=2, 3, . . . , L)) are chosen to minimize the bandwidths ofthe intermediate images, in order that they may be sampled mostsparsely. Given any such set {α_(l,m)}, the ν₂'s (shear-factors) inlevel l=1 can always be chosen to ensure that the backprojected image{circumflex over (f)} produced by the entire hierarchical algorithm ofFIG. 13 is equal to the backprojection expression (23). In particular,the shear-scale parameters in the initial level of the hierarchy dependon the set {α_(l,m)} as follows:$\sigma_{1}^{p} = \left\{ {{\begin{matrix}{\cos\quad\theta_{p}} & {{{if}\quad p} \leq \frac{P}{2}} \\{\cos\left( {\theta_{p} - \frac{\pi}{2}} \right)} & {{{if}\quad p} > \frac{P}{2}}\end{matrix}\sigma_{2}^{p}} = \left\{ \begin{matrix}{{\sin\quad\theta_{p}} - {\cos\quad\theta_{p}{\sum\limits_{l = 2}^{L}\quad{\alpha_{l,{\mu{({p,l})}}}\text{:}}}}} & {p \leq \frac{P}{2}} \\{{\sin\left( {\theta_{p} - \frac{\pi}{2}} \right)} - {{\cos\left( {\theta_{p} - \frac{\pi}{2}} \right)}{\sum\limits_{l = 2}^{L}\quad{\alpha_{l,{\mu{({p,l})}}}\text{:}}}}} & {p > \frac{P}{2}}\end{matrix} \right.} \right.$where μ(p, l)=┌p/3^(l−1)┐ describes the path in the hierarchy from thep^(th) projection in level 1 to the root.

By examination of the upper half p≦P/2) of the block-diagram in FIG. 13,following the sequence of shears backwards from the last level, it canbe seen that in this top half of the hierarchy the underlying continuousintermediate image I_(l,m) is related to its associated virtualintermediate image Ĩ_(l,m) by a shear i.e.,I _(l,m) =S _(x)(β_(l,m))Ĩ _(l,m)   (29)

Similarly to the derivation of Equation (12), it can be shown that theintermediate shear-factors depend on the {β_(l,m)} as follows:α_(1,3m−1)=β_(l+1,m)−β_(l,3m−1) for i=0, 1, 2.   (30)The freedom in the choice of the shear coefficients α_(l,m) provides afreedom in the choice of the β parameters, which can be used to minimizethe bandwidth of the intermediate projections. This will reduce theirsampling requirements, and improve the computational efficiency of thealgorithm.

Let the spectral support of the intermediate image I_(l,3m−1) be denotedby W_(l,3m−i), and its y bandwidth (i.e., slow bandwidth) be denoted byΩ_(y)(W_(l,3m−i)). The optimal β_(l,3m−i) minimizing Ω_(y)(W_(l,3m−1))is then, for i=0,1, 2 $\begin{matrix}{\beta_{l,{{3m} - i}}^{*} = {\underset{\beta \in {\mathbb{R}}}{\arg\quad\min}\left\{ {\Omega_{y}\left( {{\mathcal{S}_{y}(\beta)}{\overset{\sim}{W}}_{l,{{3m} - i}}} \right)} \right\}}} & (31)\end{matrix}$

If the view angles are uniformly distributed, it is found thatβ_(l,3m−1)*≈B_(l+1,m)*. An advantageous choice is then β_(l,3m−1)*

β_(l+1,m), yielding α_(l,3m−1)*=0, which eliminates about one third ofthe required shear operations.

The set {β_(l,m)} is determined by tracking the spectral support of theintermediate images in the algorithm, backwards from the final to firstlevel of the hierarchy. Following from the slice projection theorem, itis clear that the spectral support of Ĩ_(l,m) is{tilde over (W)} _(l,m)={(r cos θ_(p) , r sin θ_(p)):p ε N _(l,m) ,r ε[−π, π]}  (32)For the sake of sampling requirements, the spectral support {tilde over(W)}_(l,m) is characterized by its end points denoted as {tilde over(E)}_(l,m),{tilde over (E)} _(l,m)={(cos θ_(p), sin θ_(p)): p ε N _(l,m)}  (33)and shown in FIG. 12. The upper and lower band-edges of this set {tildeover (E)}_(l,m) are denoted by {tilde over (E)}_(l,m) ⁺ andrespectively. So if N_(l,m)={b, b+1, . . . , e}, {tilde over (E)}_(l,m)⁺=(cos θ_(e), sin θ_(e)) and {tilde over (E)}_(l,m) ⁻=(cos θ_(b), sinθ_(b)). Because of the ternary combination rule, the upper, middle andlower third sets of these spectral-support end-points are denoted as{tilde over (E)}_(l+1,m) ⁰, {tilde over (E)}_(l+1,m) ¹, and {tilde over(E)}_(l+1,m) ² i.e. {tilde over (E)}_(l+1,m) ^(i)={tilde over(E)}_(l,3m−i). And following from (30) and Fourier theory about affinetransformations, E_(l,m)=S_(y)(β_(l,m)){tilde over (E)}_(l,m). Theoptimal shear factors are then $\begin{matrix}{{\alpha_{l,{{3m} - i}}^{*} = {{{- \frac{\left( E_{{l + 1},m}^{i} \right)_{y}^{+} + \left( E_{{l + 1},m}^{i} \right)_{y}^{-}}{\left( E_{{l + 1},m}^{i} \right)_{x}^{+} + \left( E_{{l + 1},m}^{i} \right)_{x}^{-}}}\quad i} = 0}},2} & (34)\end{matrix}$

The algorithm FINDALPHAS shown in FIG. 14 finds the optimalshear-factors α_(l,m)* while traversing through the hierarchy down fromlevel L to level 2 in the shear-scale backprojection algorithm (FIG.13). The inputs to the algorithm are the sets of end-points of thespectral support of the images at the start of the (L+1)-th level:E_(L+1,m) for m=1, 2. In particular, in the case of uniformlydistributed angles, E_(L+1,1)=E_(L+1,2)={(cos θ_(p), sin θ_(p)): p=1, 2,. . . , P/2}.

Finding Upsampling Factors

The y-upsampling factors {U_(l,m)} in the algorithms and differentembodiments of this invention determine the sampling density of theintermediate images in the slow direction. These upsampling factors canbe chosen to meet the sampling requirements of the intermediate images,while being restricted to integer values. This restriction to integershas computational advantages, and may be expanded to that of rationalnumbers with low-valued denominators for the same reason.

The problem of choosing these factors is slightly simpler in the case ofthe shear-scale algorithm than the two-shear case. In the former case,examining FIG. 13 it is easy to see that the slow sampling interval ofthe intermediate image I_(l,m) is Π_(l′=l) ^(L) U_(l′,μ(l′,l,m)), whereμ(l′, l, m)=[m/3^(l′−l)] describes the path from I_(l,m) to the finalnode of the algorithm. Sampling theory requires that slow directionsampling interval be as follows: $\begin{matrix}{\Delta_{s}^{l,m} = {{\prod\limits_{l^{\prime} = l}^{L}\quad U_{l^{\prime},{\mu{({l^{\prime},l,m})}}}} \leq {\pi/{\Omega_{y}\left( W_{l,m} \right)}}}} & (35)\end{matrix}$where Ω_(y)(W_(l,m)) is the y-bandwidth of I_(l,m). The computationalcost of the whole algorithm, given a set of upsampling factors U

{U_(l,m) εZ:l=2, 3, . . . , L; m=1, 2, . . . , 6·3^(L−1)}, can be shownto be $\begin{matrix}{{\mathcal{J}(\mathcal{U})} = {{\sum\limits_{l = 2}^{L}\quad{\sum\limits_{m = 1}^{6 \cdot 3^{L - l}}\quad\left( \frac{c_{l,m}}{\prod\limits_{l^{\prime} = l}^{L}\quad U_{l^{\prime},{\mu{({l^{\prime},l,m})}}}} \right)}} + c}} & (36)\end{matrix}$where the constants c_(l,m) denote the relative computational cost ofdigital coordinate transformation applied to intermediate digital imageI_(l,m) ^(d). These constants are determined by the order of digitalfilters used, and by the particular implementation of the coordinatetransformation, as will be described in the discussion of efficientimplementations of coordinate transformations. The best set U tominimize J(U) subject to the constraint Equation (35), can be solved bydynamic programming or a comprehensive search. For a given set of viewangles, the best set of upsampling factors can be precomputed and storedin a lookup table.

In the case of the two-shear rotation algorithm, the problem iscomplicated slightly by the fact that the use of two-shear rotationcauses a defacto fractional up and down sampling of intermediate images.In particular, the slow-sampling periods Δ_(l,m) in adjacent levels ofthe algorithm are related as follows: Δ_(s) ^(l+1[m/3])=Δ_(s)^(l,m)/(U_(l,m) _(K) _(l,m)) where$\kappa_{l,m} = \left\{ \begin{matrix}{1\quad{if}} & {{m = 2},5,8,\ldots\quad,{{2*3^{k}} - {1\left( {{no}\quad{rotation}} \right)}}} \\{{1/{\cos\left( {\Delta_{\theta}*3^{l - 1}} \right)}}{\quad\quad}{if}} & {{m = 1},3,4,6,\ldots\quad,{2*3^{k}\left( {{{rotation}{\quad\quad}{by}} \pm {\Delta_{\theta}*3^{l - 1}}} \right)}}\end{matrix} \right.$So the constraint on the slow direction sampling interval becomes asfollows: $\begin{matrix}{\Delta_{s}^{l,m} = {{\prod\limits_{l^{\prime} = l}^{L}\quad\left( {U_{l^{\prime},{\mu{({l^{\prime},l,m})}}}\kappa_{l^{\prime},{\mu{({l^{\prime},l,m})}}}} \right)} \leq {\pi/{{\Omega_{y}\left( W_{l,m} \right)}^{\quad}}^{\quad}}}} & (37)\end{matrix}$where Ω_(y)(W_(l,m)) is the y-bandwidth of I_(l,m). The computationalcost given by Equation (36) now can be minimized subject to theconstraint in Equation (37).Oversampled Versions of the Algorithms

Oversampling can be incorporated into this algorithm to improve theaccuracy of the reconstruction. One way to do so is uniformly over allintermediate images in the algorithm: whenever an interpolation isperformed in the algorithm, the image being operated upon is at leastoversampled by some predetermined constant γ. In particular, the ratioof the Nyquist rate to the sampling frequency (in both slow and fastdirections) of every intermediate image that is subject to a digitalcoordinate transformation or resampling should preferably be less thanγ. This oversampling is preferably incorporated while modifying thealgorithms of the present invention as little as possible.

Oversampling in the Slow Direction

In the previously described algorithms of the present invention, thesampling frequency in the slow direction is controlled by the upsamplingfactors U_(l,m). This proves to be a useful tool for maintaining theoversampling condition in the oversampled versions of the algorithms,and therefore results in no alterations to the structures of thealgorithms of the present invention. The upsampling by a factor of 1/γis achieved by simply modifying the constraint on the slow-directionsampling interval by the factor γ, i.e., requiring Δ_(s)^(l,m)<γπ/Ω_(s)(I_(l,m)) for l=2, 3, . . . , L.

Oversampling in the Fast Direction

In the fast direction, the computationally inexpensive integerupsampling is not used to control the oversampling, so the algorithm ismodified to involve fractional upsampling. In the two-shear hierarchicalbackprojection algorithm, such a fractional upsampling by factors{U_(l,m): m=1, 2, . . . , 2.3^(L)} in the slow direction is alreadyincluded in level l=1. We therefore simply increase the upsamplingfactors U_(l,m) to incorporate this oversampling and then downsample theimage after the last digital coordinate transformation has beenperformed, to return the image to the desired sampling scheme (whereΔ_(f)=Δ_(s)=1.0). This modification to the two-shear algorithm thereforeinvolves only one additional level of fractional resampling. This xcoordinate resampling is combined with the digital x-shear in the L^(th)level for improved computational efficiency. The block diagram of thisalgorithm is in FIG. 15.

In the shear-scale hierarchical backprojection algorithm, thex-upsampling operations are combined with the x-shear-scales (C({rightarrow over (σ)}) at the beginning of the algorithm) to avoid an extralevel of resampling. It is easy to see that x-upsampling is simply aspecial case of x-shear-scaling, with the shear factor set to 0. Thecombination of these two operations is still an x-shear-scale: C(σ₁,σ₂)C(1/U, 0)=C(σ₁/U, σ₁/U). The downsampling at the end of the algorithmis combined with the x-shears in L^(th) level (S_(α) _(L,m)|_(m=1, . . . , 6)), to make them effectively shear-scales of thedigital image. It can be shown that an x-shear by a followed by ax-downsampling by U is effectively an x-shear-scale: C(U,0)S_(x)(α)=C(U, 1+Uα).

The exact values of these upsampling and downsampling factors isdetermined both by the parameter γ and the spectral structure of theintermediate images. In the fast direction, the oversampling conditionis that the fast direction (y coordinate) sampling interval satisfyΔ_(f) ^(l,m)<γπ/Ω_(ƒ) (I_(l,m)) for l=2, 3, . . . , L. So the additionalupsampling necessary is $\begin{matrix}{\eta = \frac{\max_{2 \leq l \leq L}{\Omega_{f}\left( I_{l,m} \right)}}{\gamma\pi}} & (38)\end{matrix}$

Consequently, the upsampling factors Ũ_(l,m), of the oversampledalgorithm, are modified from U_(l,m), those of the non-oversampledtwo-shear algorithm, as follows:Ũ _(1,m) =ηU _(1,m) and Ũ_(L+1,m)=η  (39)

In other words, the upsampling factors in the first level are modifiedto ensure that all the intermediate images in the algorithm areoversampled according to the parameter γ. Consequently, after the lastrotation, the images have a fast-sampling interval of 1/ρ; and,therefore, at the L^(th) level, they have to be downsampled by η toreturn them to a unit-sampling scheme.

It must be noted here that because the input projections are assumed tobe sampled exactly at the Nyquist rate, the oversampling condition inthe fast direction will not be satisfied for the first-level images,B₀q_(θ) _(p) , p=1, . . . , P.

The block diagram of this ternary oversampled shear-scale-basedalgorithm is in FIG. 16.

Collapsing Sequences of Similar Operations

Whenever there is a sequence of two operations operating on the samesingle coordinate, they may be combined for improved computational costand resampling accuracy. The following is another example, in additionto the previously mentioned ones of combining x-shears andx-upsampling/downsampling, or x-shear-scales withx-upsampling/downsampling.

The non-oversampled two-shear algorithm of FIG. 10 involves a y-shearfollowed by an x-shear of four out of six of the intermediate images inthe L^(th) stage. We incorporate the final downsampling of the image inthe oversampled version of FIG. 15 by collapsing the sequence of twooperations—the x-downsampling followed by the x-shear of this stage—intoa single one. This leaves the length of the cascade of interpolationsunchanged from the non-oversampling case.

Hierarchies of Arbitrary Radixes/Branching Factors

All these algorithms can be easily modified for the case when the set ofview-angles is not of the form 2*3^(L). Though the preferred embodimentis for all the branches of the hierarchy to involve the aggregation oftriplets of intermediate images, or use rotations by ±π/2 or 0,arbitrary numbers of intermediate images may be combined at each stage.

Given an arbitrary number (say M) of intermediate images at a particularlevel of the algorithm, 3×└M/3┘ of them may be combined in groups ofthree (where └x┘ is the largest integer less than or equal to x). If theremaining number of intermediate images is two then they may beaggregated as a pair to produce an image at the next level. If there isonly a single intermediate image remaining it may be passed on, withoutalteration, to the next level of the hierarchy.

The branching factor of the hierarchy (the number of branches thataggregate at a node of the hierarchy) may be altered, not just toaccommodate arbitrary numbers of view-angles, but also to reduce thedepth of the hierarchy and thereby improve image quality. In that caseit may be useful to aggregate images not in pairs or triplets but inlarger groups.

The previous prescriptions for the parameters of the coordinatetransformations may be easily extended to nodes with arbitrary branchingfactors. For example in the rotation-based algorithm whereI_(l+1,m)=Σ_(i=1) ^(M)(δ_(l,m) _(i) )I_(l,m) _(i) , the relationsdescribed in Equations (10) and (15) still hold, and the rotation anglestherefore are prescribed by δ_(l,m) _(i) =α_(l+1,m)−α_(l,m) _(i) . Inthe case of the shear-scale based algorithm where I_(l+1,m)=Σ_(i−1) ^(M)S_(x)(α_(l,m) _(i) )I_(l,m) _(i) , Equation (29) still holds. Extendingthe notation E_(l,m) ^(i) to refer to the i^(th) set of the M sets ofprojections being aggregated, one obtains an equation for α_(l,m) _(i) *identical to the right-hand-side of Equation (34).

0.1 Hierarchical Algorithms Based on Other Image Transformations

The back-projection equation (3) may be written in the form$\begin{matrix}{{\hat{f}\left( \overset{\rightarrow}{x} \right)} = {\frac{1}{P}{\sum\limits_{p = 1}^{P}{\left\lbrack {\mathcal{B}_{0}q_{p}} \right\rbrack\left( {A_{\theta_{p}}\overset{\rightarrow}{x}} \right)}}}} & (40)\end{matrix}$for any matrix A_(θ), as long as the first row of A_(θ) is [cos θ sinθ]. The freedom to choose arbitrary values for the remaining two entriesof A_(θ) allows for flexibility in the design of the coordinatetransformations used in the hierarchical algorithm. Matrix A_(θ) can befactored as A_(θ)=Π_(l=1) ^(L) A(δ_(l)) for some parameter vectors δ_(l)that are related to θ, but are not completely determined, owing to thefreedom in the two bottom entries of A_(θ). This factorization can beused to derive a hierarchical decomposition of the backprojectionequation (40) with a corresponding block diagram such as FIG. 5, withthe coordinate transformation steps denoted by K(δ_(l,m)) representingthe image coordinate transformation operators defined as(K(δ_(l,m))f)({right arrow over (f)})=ƒ( A(δ_(l,m)){right arrow over(x)}).

Clearly, the specific embodiments of this invention described herein arespecial cases of this more general choice of matrix A(o and itsfactorization, and the associate coordinate transformation. The effectof these coordinate transformations on the Fourier spectrum of theintermediate images is analyzed similarly to the cases alreadydescribed, because the effect of an affine transformation by matrix A inthe spatial domain is an affine transformation by matrix A^(−T), i.e.,the inverse-transpose of A, in the frequency domain. Similarconsiderations can therefore be used to select free parameters in thetransformations with the goal of reducing the computationalrequirements. Thus, the class of digital image coordinatetransformations used in the hierarchical backprojection algorithms ofthe present invention includes many other transformations in addition tothose described for specific embodiments. Furthermore, because matricesA(δ_(l,m)) can be factorized into a product of triangular matrices, thecoordinate transformations can be performed as a cascade of single-axiscoordinate transformations, if desired.

0.2 Efficient and Accurate Implementation of Digital Image CoordinateTransformations and Resampling

The accuracy and speed of the hierarchical backprojection andreprojection algorithms of the present invention depend on the specificsof the implementation of the various digital image coordinatetransformations and resampling operations. Improved accuracy requires ingeneral high order filters or interpolations, which usually increasesthe computation cost. The cost of high order filtering or interpolationcan be reduced by decomposing the operations into lower orderoperations. Additional reduction in computation and/or memoryrequirements is obtained if the filters used are shift invariant. Highorder finite impulse response filters can be implemented efficientlyusing low-order recursive filter, or by using the fast Fourier transform(FE”).

In particular, if a separable representation basis is assumed for thecontinuous image, a digital y-shear of the digital image can be achievedby fractional delays of the individual vertical lines in a digital imagearray. Likewise, a digital x-shear can be expressed as a fractionaldelay of one row at a time. In turn, a fractional delay of a 1D signalcan be accomplished using a shift-invariant filter. Similardecompositions are known for 3D images and shear operations.

Resampling a digital image can also be usually performed using lowerdimensional operations, which can often be shift invariant. For example,image upsampling along one coordinate by an integer factor U may bedecomposed into U different, computationally efficient, fractionaldelays. This is essentially, the so-called polyphase decomposition,well-known in digital signal processing. Rational, but non-integerresampling along one coordinate can be decomposed into a cascade ofinteger down and upsampling, each of which is efficiently performed.

More general digital image resampling can also be decomposed into lowerdimensional operations. Consider the resampling of a digital 2D image ƒfrom a sampling pattern with samples S₁ lying on a family of curves,denoted CF₁, to another another pattern with points S₂ lying on adifferent family of curves, CF₂, producing the digital image h. If thetwo families of curves intersect at sufficient density the method thatwas described with reference to FIG. 23 for resampling from one fan toanother rotated fan can be used for general curves. Otherwise a thirdfamily of curves CF₃ can be introduced, which intersects both CF₁ andCF₂ at a desired density. The digital image can then be resampled fromCF, to its intersections with CF₃, then to the intersections of CF₃ withCF₂, and finaly on CF₂ to the desired sampling pattern.

This process generalizes to 3D, for example by considering surfacesinstead of curves, to first reduce the process to one of resampling onsurfaces, and then using the 2D process to further reduce resampling onsurfaces to resampling on curves.

Digital coordinate transformations and resampling can often be combined,improving the computational efficiency. For example, the digitalx-shear-scale operations used in the shear-scale algorithm shown in FIG.13 can be decomposed into resampling operations on individual horizontallines.

1 Divergent-Beam Fast Hierarchical Backprojection Algorithms

1.1 Fan-Beam Projection and Backprojection

Consider the case of equiangular-spaced detectors located on a circlearound the object. The fan-beam tomographic projection, at asource-angle β, of a two-dimensional image f(x, y) is denoted by(Rβf)(γ) and is defined as the set of line integrals along the rays of afan, parameterized by γ, centered at the source position on a circle ofradius D from the origin. The function ƒ( {right arrow over (x)}) isassumed to be zero-valued outside a disc of radius R.

The fan-beam ray V(β, γ, T), shown in FIG. 17, is defined by$\begin{matrix}{{V\left( {\beta,\gamma,T} \right)} = {{D\begin{bmatrix}{{- \sin}\quad\beta} \\{\cos\quad\beta}\end{bmatrix}} + {T\begin{bmatrix}{\sin\quad\left( {\beta + \gamma} \right)} \\{- {\cos\left( {\beta + \gamma} \right)}}\end{bmatrix}}}} & (41)\end{matrix}$The fan-beam projection at source-angle β and fan-angle γ is(R_(β)f)(γ)=∫_(T=0) ^(∞)ƒ( V(β,γ, T))dT. Since ƒ is zero outside thedisc of radius R, the integral needs only be performed within the disci.e. between T_(ST) and T_(END).

In computed tomography with the fan-beam geometry, projections areavailable at a set of discrete source angles {β_(p): p=1, 2, . . . , P},and within each fan the angles of the rays are indexed by {γ_(j): j=1,2, . . . , J}. In the case of equiangular fan-beam geometry thedetectors are equally distributed on the arc of a circle centered at thesource, so the fan-angles are evenly spaced. In the case of equispaceddetectors, the detectors are equally distributed on a line perpendicularto the line from the source position to the origin. The fastbackprojection algorithms for fan-beam geometry described here assumeequiangular distribution with J odd and γ_(j)=γ_(γ)·(j−(J+1)/2).However, the algorithms may be easily extended to other fan-anglegeometries.

The reconstruction algorithm from a set of P fan-beam projections {

_(β) _(i) (γ): i=1, 2, . . . , P}, may be expressed as the scaling andfiltering of each fan-beam projection followed by a weightedbackprojection (5): $\begin{matrix}{{\hat{f}\left( \overset{\rightarrow}{x} \right)} = {\sum\limits_{p = 1}^{P}{{W\left( {T\left( {\overset{\rightarrow}{x},\beta_{p}} \right)} \right)}{q_{\beta_{p}}\left( {\gamma\left( {\overset{\rightarrow}{x},\beta_{p}} \right)} \right)}}}} & (42)\end{matrix}$where W(T) is an appropriate weight function, T({right arrow over(()}x), β)) and y({right arrow over (()}x), β) are the distance alongthe ray between source and image point {right arrow over (x)} for sourceposition β, and the fan angle of that ray, respectively, and q_(β)(γ)are the weighted and filtered projectionsQ _(β)(γ)=′ _(β)(γ)*g(γ), ′ _(β)(γ)= _(β)(γ)*D*cos γwhere g(γ) is an appropriate filter.

Definition 1.1 The weighted backprojection of a function q at a singlesource angle fi is defined by(W_(β)q)({right arrow over (x)})

W(T({right arrow over (x)}, β))q(γ({right arrow over (x)}, β)) tm (43)and the weighted backprojection operator W_({right arrow over (β)})(where {right arrow over (β)} ε [0, 2π)^(P)) is defined by$\begin{matrix}{{\left( {\mathcal{W}_{\overset{\rightarrow}{\beta}}\left\{ q_{p} \right\}_{p = 1}^{P}} \right)\left( \overset{\rightarrow}{x} \right)}\overset{\Delta}{=}{\sum\limits_{p = 1}^{P}{\left( {\mathcal{W}_{\beta_{p}}q_{p}} \right)\left( \overset{\rightarrow}{x} \right)}}} & (44)\end{matrix}$Therefore {circumflex over (f)}({right arrow over(x)})=(W_({right arrow over (β)}){q_(p)}_(p=1) ^(P))({right arrow over(x)}).

It is easily shown thatW _(β) =K(−β)W ₀   (45)where K(β) denotes the rotation by β operator, and consequently that$\begin{matrix}{{\left( {\mathcal{W}_{\overset{\rightarrow}{\beta}}\left\{ q_{p} \right\}_{p = 1}^{P}} \right)\left( \overset{\rightarrow}{x} \right)} = {\sum\limits_{p = 1}^{P}{\left( {{\mathcal{K}\left( {- \beta_{p}} \right)}\mathcal{W}_{0}q_{p}} \right)\left( \overset{\rightarrow}{x} \right)}}} & (46)\end{matrix}$

Thus, as in the parallel-beam case, the weighted backprojection of Pfan-beam projections may be expressed as the sum ofweighted-zero-backprojected images as seen in FIG. 4, with B₀=W₀. Infact, close analogs of the backprojection operator defined in Equations(43) and (44) apply more generally to the reconstruction of functions intwo and higher dimensions from projections of a general form, withappropriate definition of the functions γ({right arrow over (x)}, β) andT({right arrow over (x)}, β). The methods of the present inventionextend to these other applications, with appropriate definition of acoordinate transformation K.

1.2 Fast Hierarchical Backprojection

The Fast Hierarchical Backprojection Algorithm for the fan-beam geometryis similar to that for the parallel-beam case. It combines the fan-beamprojections in a ternary hierarchical structure, exploiting the factthat intermediate images formed by projections that are close to eachother in projection angle fi can be sampled sparsely. The block-diagramis shown in FIG. 18. Similar algorithms can be derived for binary orother (possibly mixed) radix structures.

The equations that govern the combination of the underlying continuousimages in FIG. 18 are as follows:I_(1,m)

W₀q_(m)   (47)I_(l+1,m) K(δ _(l,3m−2))I_(l,3m−2) +I _(l,3m−1) +K(δ_(l,3m))I_(l,3m)l=1, 2, . . . , L   (48)I_(L+1,1)

I_(L,1)+K(−π/2)I_(L,2)+K(−π)I_(L,3)+K(−3π/2)I_(L,4)   (49)

In the embodiment described next, it is assumed that the source anglesare uniformly spaced in [0, 2π) as follows:$\beta_{j} = {{{- \frac{\pi}{4}}\left( {1 - \frac{1}{9}} \right)} + {\Delta_{\beta}\left( {j - 1} \right)}}$where Δ_(β) = π/(2 ⋅ 3^(L))

The intermediate rotation angles are then chosen as in the parallel-beamcase using Equation (22), with Δ_(β) replacing Δ_(θ).

The constituent fans of an intermediate image in three levels of thealgorithm with block diagram FIG. 18 are illustrated in FIG. 20.

1.2.1 The Fan-Beam Sampling Scheme

One embodiment of the fan-beam algorithm uses a sampling scheme of theintermediate images in the algorithm derived by analogy to theparallel-beam case, which is described next. Alternatives andimprovements are described later.

A single-fanbeam-backprojection image (an image W₀q produced by weightedbackprojection of a single fanbeam) has a structure amenable to sparsesampling. It follows from Equation (43)that (W_(β)q) (V(β, γ,T))=W(T)q(γ). A single sample at a particular T=T₁ on the ray indexed byγ in a fan oriented at β=0 is therefore sufficient to specify the valueof the image along the entire ray: it is simply proportional to W(T)where T is the distance from the source to the point in the ray.

An intermediate image in the algorithm is a sum of several rotatatedversions of such single-fanbeam-backprojection images, each generatedfrom a constituent projection. We refer to the fan at the angle that isat the center of the angular interval spanned by the constituentprojections as the central constituent fan. This may correspond to anactual projection—e.g., in the case of a ternary algorithm,—or a tovirtual projection—e.g., in the case of a binary algorithm. For example,the fan in FIG. 20(a) is the central constituent fan for both FIG. 20(b)and FIG. 20(c).

Recall that in the parallel-beam algorithm the intermediate images aresampled on cartesian sampling patterns aligned with the centralconstituent. In the rotation-based parallel-beam case the intermediateimages are sampled along parallel vertical rays. The spacing of samplesalong each of these parallel rays is chosen to guarantee that the otherconstituent projections of the image are sufficiently sampled (by theNyquist sampling criterion) along that ray. The necessary spacing isexactly equal to the spacing of the intersection points of the verticalrays with a set of rotated parallel rays corresponding to the extremalconstituent projection, that is, the one farthest away in view anglefrom the center projection.

In direct analogy to the parallel-beam case, here the intermediateimages are sampled along the rays of the central constituent fan. Thesampling points along the rays of the central fan are the intersectionpoints of the central fan with the extremal constituent fans. Thisresults in samples that are not uniformly spaced along each ray.Consider two fans V(β, γ_(j), •)|_(j=1) ^(J) and V(β′, γ′_(j), •)|_(j=1)^(J). Assuming the equiangular fanbeam geometry with an odd number J ofdetectors, γ_(j)=Δ_(γ·(j−(J+)1)/2). The points on the r^(th) ray of theβ-fan (corresponding to fan-angle γ_(r)) that intersect the β′-fan, aredetermined by the equation V (β, γ_(r), T)=V(β′, γ′_(j), T′) which hasthe solution${T_{\beta,\beta^{\prime},\gamma_{r}}\left( \gamma_{j}^{\prime} \right)} = \frac{D\left\lbrack {{\sin\left( {\beta - \beta^{\prime} - \gamma_{j}^{\prime}} \right)} + {\sin\left( \gamma_{j}^{\prime} \right)}} \right\rbrack}{\sin\left( {\beta - \beta^{\prime} + \gamma_{r} - \gamma_{j}^{\prime}} \right)}$

The function T_(β,β′,γ) _(r) (γ′) describes how the fan V(β′, •, •)varies along the r^(th) ray of the fan V(β, •, •). It carriesinformation about the varying sampling rate along a ray of the fan. Thesteeper the slope, the sparser are the samples. The local sampling rateat any γ′ is proportional to (dT/dγ′)⁻¹. A typical T(γ′) is displayed bythe dashed line in FIG. 22.

Two modifications are made to T_(β,β′,γ) _(r) (γ′) to get a new samplingfunction {tilde over (T)}_(β,β′,γ) _(r) (γ′):

1. Intermediate images are to be sampled within the disc of radius R,with margins of at least one sample on each ray outside the disc (onesample with T<T_(ST) and one sample with T>T_(END)). Defineγ′_(END)=T_(β,β′,γ) _(r) ⁻¹ (T_(END)): the value of γ′ corresponding toT_(END). Fans with adjacent source angles may lack intersection pointswith T>T_(END). In order to rectify that, the slope of {tilde over(T)}_(β,β′,γ) _(r) (γ′) for T>T_(END) is clamped to that at T_(β,β′,γ)_(r) (γ′_(END)).

2. The final image is sampled on a Cartesian pattern with unit samplingintervals. This is the smallest interval at which intermediate imagesneed to be sampled. The point γ′* at which a unit sampling rate alongthe γ^(th) ray is achieved is determined by solving the equationd/(dγ)T_(β,β′,γ) _(r) (γ)=1/Δ_(γ) for γ. To maintain the unit samplinginterval constraint, the slope of {tilde over (T)}_(β,β′,γ) _(r) (γ′)for γ′≧γ′* is set to 1/Δ_(y).

FIG. 22 displays T_(β,β′,γ) _(r) (γ) and {tilde over (T)}_(β,β′,γ) _(r)(γ) for a typical β and γ. It also displays the T-values of the actualsampling points on the ray which are the set {{tilde over (T)}_(β,β′,γ)_(r) (jΔ_(γ)): j ε Z}. The integers j are chosen so that there aremargins of at least one sample on each side of the image boundary.

Clearly, the locations of the sample points prescribed by the principlesoutlined herein can be computed from the geometry of the fans, Equation(21) for the set N_(l,m) of the constituent projections, the selectedrotation angles δ_(l,m) for the intermediate images (given in Equation(22) for the case of equispaced view angles), and the chosen form for{tilde over (T)}_(β,β′,γ) _(r) (γ).

Separable Rotation and Up-Interpolation for Fan-Beam Sampling

As described in the Overview, the upsampling operations, and upsamplingcombined with rotation operations can be decomposed into computationallyefficient one dimensional resampling operations.

1.2.2 Sampling Scheme Based on Local Fourier Structure

Given an image ƒ( {right arrow over (x)}):

^(n)→R, we want to find a sampling function t({right arrow over (x)}):

^(n)−

^(n) such that ƒ( t({right arrow over (x)})) has a small essentialbandwidth and therefore can be sampled on the set of points {t({rightarrow over (m)}): {right arrow over (m)} ε

^(n)}. We find this sampling function as follows. As will be explainedshortly, knowing how ƒ( {right arrow over (x)}) is composed of itsconstituent projections, we find the matrix function${\upsilon\left( \overset{\rightarrow}{x} \right)} = \begin{bmatrix}{\upsilon_{11}\left( \overset{\rightarrow}{x} \right)} & {\upsilon_{12}\left( \overset{\rightarrow}{x} \right)} \\{\upsilon_{21}\left( \overset{\rightarrow}{x} \right)} & {{\upsilon_{22}\left( \overset{\rightarrow}{x} \right)}\quad}\end{bmatrix}$that describes how the function ƒ should be sampled locally at the point{right arrow over (x)}. We then integrate this matrix function over theimage domain to get the sampling function t({right arrow over (x)}).

In our algorithms, we know that the intermediate image ƒ is of the form:f=Σ_(p)K(δ_(p))W₀q_(p) for some angles δ_(p) ε [β_(min), β_(max)].Ignoring the weighting by W(T), and assuming that the projections q_(p)are sampled exactly at the Nyquist rate, we know the local structure ofthe spectral support at each point in the image ƒ. Because of thefan-backprojection of each band-limited projection q_(p), an image thatis fan-backprojected from a single projection has a negligible spatialbandwidth in the direction of the spoke of the fan, while it's bandwidthin the perpendicular direction is inversely proportionate to thedistance from the vertex of the fan as seen in FIG. 24. When a set ofthese fans are rotated and added together, the local spectral support ofthis resulting image is the union of the spectral supports of theindividual fans.

For example when the constituent fans have source angles between β_(min)and β_(max), as shown in FIG. 25, the local spectral support at a point{right arrow over (x)} in the image domain is a warped wedge with radiioriented between θ_(min) and θ_(max), and a radial bandwidth of ΩR/r_(θ)at angle θ. Here Ω is the spatial bandwidth of the back-projectedprojection at the center of the image plane assuming the projections aresampled at the Nyquist rate. In the equiangular case Ω is number ofsamples of the projection divided by the length of the arc through theorigin of the image plane that the backprojected projection covers.Mathematically the spectral support is $\begin{matrix}{\left\{ {{\left( {{\cos\quad\theta\quad\omega_{\theta}},{\sin\quad\theta\quad\omega_{\theta}}} \right)\text{:}\theta}\quad \in {\left\lbrack {\theta_{\min},\theta_{\max}} \right\rbrack\quad{and}\quad\omega_{\theta}} \in \left\lbrack {{- \Omega_{\theta}},\Omega_{\theta}} \right\rbrack} \right\}{where}{\theta_{\max} = {\tan\frac{{R\quad\sin\quad\beta_{\max}} - x_{1}}{{R\quad\cos\quad\beta_{\max}} + x_{2}}}}{\theta_{\min} = {\tan\frac{{R\quad\sin\quad\beta_{\min}} - x_{1}}{{R\quad\cos\quad\beta_{\min}} + x_{2}}}}{\Omega_{\theta} = {R\quad{\Omega/r_{\theta}}}}{r_{\theta} = \sqrt{{\left( {x_{1} - {R\quad\sin\quad\theta}} \right)^{2} + \left( {x_{2} - {R\quad\cos\quad\theta}} \right)^{2}}\quad}}} & (50)\end{matrix}$

Given this knowledge of the local two-dimensional fourier structure at apoint {right arrow over (x)} in the image, the matrix function v({rightarrow over (x)}) at that point is the sampling matrix that, if thespectral support at that point were uniform across the whole image (suchas in the parallel-beam case), would efficiently sample the image. Inthe intermediate images of interest here, this produces two distinctsmall-bandwidth and large-bandwidth sampling directions. We fix thefirst coordinate to be the small-bandwidth direction.

We approximate this sampling matrix function v({right arrow over (x)})by $\left( {{\begin{bmatrix}s_{1} & 0 \\0 & s_{2}\end{bmatrix}\begin{bmatrix}1 & \alpha_{1} \\0 & 1\end{bmatrix}}\begin{bmatrix}1 & 0 \\\alpha_{2} & 1\end{bmatrix}} \right)^{T}.$Here $\begin{matrix}{{\alpha_{2} = {{- \frac{{\Omega_{\min}\sin\quad\theta_{\min}} + {\Omega_{\max}\sin\quad\theta_{\max}}}{{\Omega_{\min}\cos\quad\theta_{\min}} + {\Omega_{\max}\cos\quad\theta_{\max}}}}\quad{and}}}{\alpha_{1} = \frac{{\Omega_{\min}\cos\quad\theta_{\min}} + {\Omega_{\max}\cos\quad\theta_{\max}}}{{\Omega_{\max}\sin\quad\theta_{\min}} + {\Omega_{\min}\sin\quad\theta_{\min}}}}{s_{1} = {\left( {{{{\left\lbrack {1\quad 0} \right\rbrack\begin{bmatrix}1 & \alpha_{1} \\0 & 1\end{bmatrix}}\begin{bmatrix}1 & 0 \\\alpha_{2} & 1\end{bmatrix}}\begin{bmatrix}{\cos\quad\theta_{mid}} \\{\sin\quad\theta_{mid}}\end{bmatrix}}\Omega_{mid}} \right)^{- 1}\quad{and}}}{s_{2} = \left( {{{\left\lbrack {0\quad 1} \right\rbrack\begin{bmatrix}1 & 0 \\\alpha_{2} & 1\end{bmatrix}}\begin{bmatrix}{\cos\quad\theta_{\max}} \\{\sin\quad\theta_{\max}}\end{bmatrix}}\Omega_{\max}} \right)^{- 1}}} & (51)\end{matrix}$

Here Ω_(min)

Ω_(θ) _(min) and θ_(max)

Ω_(θ) _(max) . Angle θ_(mid) refers to the angle of the projectioncorresponding to the source-angle (β_(min)+β_(max))/2, and Ω_(mid)

Ω_(θ) _(mid) . These parameters, chosen by geometrical arguments on thespectral support, ensure that the spectral support upon transformationby the sampling matrix is restricted to [−π, π]×[−π, π].

The sampling function for the whole image is found by integrating thissampling matrix function across the whole image—solving the set ofdifferential equations: dt_(i)/dn_(j)=v_(ij)(t({right arrow over (n)}))for i, j=1, 2. This may be solved numerically. Even if this exactprescribed pattern is not used, the local Fourier support analysis willevaluate the effectiveness of any sampling pattern.

The resulting sampling patterns for a few intermediate images are shownin FIGS. 26(A) and 26(B). For the fan-beam case, this localFourier-based method produces sampling patterns that are similar to theones resulting from the intersection-based methods described earlier.dThe previously described separable method to resample from one samplingpattern to another can be used with these patterns also. This localFourier sampling method may be applied directly to find sampling schemesfor arbitrary projection geometries over lines, curves or planes overarbitrary dimensions. In the case of the parallel-beam geometry, itreduces to that discussed earlier.

1.2.3 Alternative Sampling Schemes

The sampling schemes in which the samples are located on a fan whosevertex is on the source trajectory (called a “sampling fan”) has someshortcomings. The sampling points are chosen separately for each rayresulting in a scheme that is not necessarily optimal in a twodimensional sense. Though the final image in the algorithm is sampled ona uniform rectangular grid with unit-spacing, the intermediate imagesusing the above mentioned scheme are sampled more densely than necessaryin certain regions of the image (eg. close to the vertex of the fan). Inorder to rectify this the above described scheme might be modified tomake sparse the sampling in such oversampled regions. Two suchpossibilities are illustrated in FIG. 27. Both these possibiliesincorporate the fact while in the first level the image is sampledefficiently on a fan, in the final level it needs to be sampled on arectangular grid. These schemes attempt to incorporate the gradualtransition from sampling on a fan to a grid.

In FIG. 27(A) samples are selected on a pseudo-fan whose vertex islocated further from the origin than the source radius. In successivelevels intermediate images are formed from fans from a larger range ofsource angles and the distance of the vertex of the pseudo-fan from theorigin increases. In the final level the vertex is at infinity; i.e. therays are the parallel lines of the rectangular grid. In FIG. 27(B) thesamples are located, not on a fan, but on a beam that is less divergentnearer the source position. In successive levels the beam becomes moreparallel and less divergent. These or other sampling schemes that take atwo dimensional (or, for 3D images, a 3D) point of view will contributeto a faster algorithm.

1.2.4 Other Optimizations

In the fan-beam case with a full trajectory, source angles aredistributed in [0, 2π). This allows us to take advantage of rotations by−π/2, −π and −3π/2 that involve no interpolation but only a rearrangingof pixels.

* The positions of points in the sampling patterns for all levels of thealgorithm can be pre-computed and stored in a look up table, or computedon-the-fly along with the processing of data. In either case,computation and/or storage can be reduced by taking advantage of thesmooth variation of sample position as a function of ray index, whichmay be observed in FIGS. 23(A) and 23(A).

1.2.5 Oversampled Fast Hierarchical Backprojection

As in the parallel-beam case, oversampling by a factor γ<1.0 can be usedto improve accuracy. One way to achieve this in the fan-beam case, is todetermine the denser sampling patterns using a fan-angle spacingΔ′_(β)=Δ_(β)·γ.

1.2.6 Algorithms for Short Scan and Alternative Angle Distributions

The methods of this invention are also applicable to the so-calledshort-scan data acquisition format, where the range of source angles issmaller than [0, 2π). The difference is in the weighting used inbackprojection: an extra weighting factor (e.g., Parker weighting) isintroduced to account for the short scan. Likewise, the algorithmsgeneralize to nonuniform distributions of the source angle, or evensmall sets of source angles that may not allow reconstruction of acomplete image. Furthermore, the order of operations in the algorithmmay be modified in various ways (for example pipelining), which would beknown to those skilled in the art. For example, the operations may beordered such that processing may begin as soon as a single projection isavailable.

1. A method (FIG. 5A) for creating a pixel image {circumflex over (f)}from projections (q₁ . . . q_(p)) comprising the steps of: (a) producing(100) intermediate-images (I₁, P) from selected projections (q₁ . . .q_(p)); (b) performing digital image coordinate transformations (102) onselected intermediate-images (I₁, P), the parameters of coordinatetransformations being chosen to account for view-angles of theprojections from which the intermediate images have been produced, andfor the Fourier characteristics of the intermediate-images; (c)aggregating subsets of the transformed intermediate-images (104)produced in step (b) to produce aggregate intermediate-images (I₂, P/2);and (d) repeating steps (b), and (c) in a recursive manner until all ofthe projections and intermediate images have been processed andaggregated to form the pixel image {circumflex over (f)}; wherein thecoordinate transformation parameters are chosen so that the aggregatesof the intermediate-images (104) may be represented with desirableaccuracy by sparse samples.
 2. The method in claim 1, in which saidaggregation (104, 108) is performed by adding digital images.
 3. Themethod in claim 1, wherein at least some intermediate images (I_(n), m)are produced in step (a) by backprojecting selected projections (q₁ . .. q_(p)).
 4. The method of claim 1 wherein at least some intermediateimages (I_(n), m) are each formed by backprojecting two or more selectedprojections (q₁ . . . q_(p)) in step (a).
 5. The method of claim 1wherein at least some aggregate intermediate images (I_(n), m) are eachformed by aggregating three or more selected transformed intermediateimages in step (d).
 6. The method of claim 1, wherein the digital imagecoordinate transformations are performed using digital filtering.
 7. Themethod of claim 1, wherein selected coordinate transformations includedigital image rotations.
 8. The method of claim 1, wherein selectedcoordinate transformations include digital image shearing (FIG. 10B,120, 122), or shear-scaling.
 9. The method (FIG. 15) of claim 1, whereinselected coordinate transformations include upsampling (101, 106) and/ordownsampling (109) of the digital images.
 10. The method of claim 8 inwhich said digital image shearing is performed by one-dimensional lineardigital filters.
 11. The method of claim 9 in which said digital imageupsampling and/or downsampling are performed by one-dimensional lineardigital filters.
 12. The method of claim 10 in which at least some ofsaid digital filters are shift-invariant.
 13. The method in claim 11 inwhich at least some of said digital filters are shift-invariant.
 14. Themethods in claims 12 or 13, in which at least some of said digitalfilters are recursive.
 15. The methods in claims 12 or 13, in which atleast some of said digital filters are implemented using a fast Fouriertransform (FFT).
 16. The method in claim 1, in which selectedoversampling is applied to selected intermediate images and/ortransformed intermediate images and/or aggregate intermediate images.17. The method in claim 1, in which non-Cartesian sampling patterns areused.
 18. The method in claim 1, in which selected coordinatetransformations may be combined within a level, or across adjacentlevels of the hierarchy.
 19. A method (FIG. 31) for creating a pixelimage {circumflex over (f)} from projections (q₁ . . . q_(p)) along acollection of lines, curves, or surfaces comprising the steps of: (a)producing (184) intermediate images (I_(l), m); (b) performing digitalimage resampling on selected intermediate images (186), the location ofsamples being chosen to account for the view-angles of the selectedprojections and for the Fourier characteristics of the intermediateimages, (c) aggregating (190) selected subsets of the resampledintermediate-images to produce aggregate intermediate-images (I_(x), m);and (d) repeating steps (b) and (c) in a recursive manner, at each levelof the recursion increasing the density of samples of the intermediateimages, until all of the projections and intermediate images have beenprocessed and aggregated to form the pixel image; wherein the samplingscheme is chosen so that aggregates of the resampled intermediate-imagesmay be represented with desirable accuracy by sparse samples.
 20. Themethod of claim 19, wherein at least some intermediate images areproduced in step (a) by weighted backprojection (FIG. 19, 180, 182) ofselected projections.
 21. The method of claim 19 wherein at least someintermediate images are each formed by weighted backprojection (FIG. 19,180, 182) of two or more selected projections in step (a).
 22. Themethod of claim 19 wherein at least some aggregate intermediate imagesare each formed by aggregating three or more selected transformedintermediate images in step (d).
 23. The method of claim 19, in whichthe intermediate images have samples that lie on a family of lines,curves or surfaces.
 24. The method of claim 19, in which the digitalimage resampling is performed by a sequence of lower-dimensional digitalfiltering operations by utilizing intermediate sampling schemes that lieon the intersections of the families of lines, curves or planes.
 25. Themethod of claim 19, in which a selected degree of oversampling isapplied to the selected resampled intermediate images, and aggregatedintermediate images.
 26. The method in claim 19, in which saidaggregation is performed by adding digital images.
 27. The method inclaim 19, in which the resampling and aggregation may be combined acrosssuccessive levels.
 28. The method of claim 30, in which at least oneintermediate image is weighted before and after resampling step (b) areincluded.
 29. The method of claim 19, in which changes in samplingdensity are accomplished by digital filtering.
 30. A method (FIG. 18)for creating a pixel image {circumflex over (f)} from projections (q₁ .. . q_(p)) comprising the steps of: (a) producing (99) a plurality ofintermediate-images (I_(l), m), with at least one corresponding to anon-Cartesian and/or non-periodic sampling pattern; (b) performingdigital image upsampling or downsampling (106) on selectedintermediate-images; (c) performing digital image coordinatetransformations on upsampled/downsampled intermediate-images; (d)aggregating (110) subsets of the transformed intermediate-imagesproduced in step (c) to produce aggregate intermediate-images; and (e)repeating steps (b), (c) and (d) in a recursive manner until all of theprojections and intermediate images have been processed and aggregatedto form the pixel image; wherein at least one of the digital imagecoordinate transformations is performed with a non-Cartesian and/ornon-periodic sampling pattern, and the coordinate transformationparameters are chosen so that the aggregates of the intermediate-imagesmay be represented with desirable accuracy by sparse samples.
 31. Themethod in claim 30, in which said aggregation is performed by addingdigital images.
 32. The method in claim 30, wherein at least someintermediate images are produced in step (a) by backprojecting selectedprojections.
 33. The method of claim 30 wherein at least someintermediate images are each formed by backprojecting two or moreselected projections in step (a).
 34. The method of claim 30 wherein atleast some aggregate intermediate images are each formed by aggregatingthree or more selected transformed intermediate images in step (c). 35.The method of claim 30, wherein the digital image coordinatetransformations are performed using digital filtering.
 36. The method ofclaim 30, wherein selected coordinate transformations include digitalimage rotations.
 37. The method of claim 30 in which said digital imageresampling upsampling and/or downsampling are performed byone-dimensional linear digital filters.
 38. The method of claims 37 inwhich at least some of said digital filters are shift-invariant.
 39. Themethod in claim 37, in which at least some of said digital filters arerecursive.
 40. The method in claim 38, in which at least some of saiddigital filters are implemented using a fast Fourier transform (FFT).41. The method in claim 30, in which selected oversampling is applied toselected intermediate images and/or transformed intermediate imagesand/or aggregate intermediate images.